The CDF from the characteristic function

[EngWiPy]

This thread will be a collection of multiple questions I asked before over different forums. I will start from the beginning, and I hope someone will follow the steps with me, because I did it before alone, and I ended with a numerical integration that is not finite, which doesn't make sense, since the result should be between 0 and 1.

My main problem is to find the cumulative distribution function (CDF) of the following random variable (RV)

$$Y=\sum_{k=1}^K\underbrace{\frac{a_k}{b_k}}_{X_k}$$

where $\{a_k, b_k\}$ are independent and identically distributed (i.i.d.) exponential RVs with parameter 1.

At first, I wanted to use the moment generating function (MGF), but some pointed out that the RVs $X_k$ have undefined means, and the MGF doesn't work in this case. They suggested the characteristic function (CF), because the CF always exists even if the means don't.

My fist step was then to find the CF of the RV $Y$ as follows

$$\psi_Y(t)=\text{E}_Y\{e^{jtx}\}=\text{E}_{X_1,\,X_2,\ldots,\,X_K}\{e^{jt\sum_{k=1}^KX_k}\}=\prod_{k=1}^K\text{E}_{X_k}\{e^{jtx_k}\}=\prod_{k=1}^K\psi_{X_k}(t)$$

So, the next step would be to find $\psi_{X_k}(t)$, which is given by

$$\psi_{X_k}(t) = \text{E}_{X_k}\{e^{jtx_k}\} = \int_0^{\infty}e^{jtx_k}\,f_{X_k}(x_k)\,dx_k$$

where $f_{X_k}(x_k)$ is the probability distribution function (PDF) of the RV $X_k$. It can be shown that $f_{X_k}(x_k)=1/(1+x_k)^2$, and thus

$$\psi_{X_k}(t) = \int_0^{\infty}\frac{e^{jtx_k}}{(1+x_k)^2}\,dx_k$$

My first question is how to evaluate this integral?

My approach: by integration by parts, this results in

$$\psi_{X_k}(t) = 1+jt\int_0^{\infty}\frac{e^{jtx_k}}{1+x_k}\,dx_k=1+jt\int_0^{\infty}\frac{\cos(tx_k)}{1+x_k}\,dx_k-t\int_0^{\infty}\frac{\sin(tx_k)}{1+x_k}\,dx_k$$

which can be written in terms of cosine and sine integrals as

$$\psi_{X_k}(t) = 1+jt\left[-\sin(t)\text{si}(t)-\cos(t)\text{ci}(t)\right]-t\left[\sin(t)\text{ci}(t)-\cos(t)\text{si}(t)\right]$$

where $\text{si}(t)=-\int_t^{\infty}\sin(z)/z\,dz$ and $\text{ci}(t)=-\int_t^{\infty}\cos(z)/z\,dz$.

Is what I did so far correct? I appreciate if some comment on this, because apparently I made a mistake in my mathematical analysis that resulted in an indefinite integral, but not sure where my mistake is.

Thanks in advance

 

[mathman]

What indefinite integrals?

 

[EngWiPy]

I meant the integral gave undefined results. Anyway, after some mathematical manipulation, I reached to the formula

$$\psi_{X_k}(t)=1+je^{-jt}E_1(-jt)$$

where $E_1(x)=\int_x^{\infty}\frac{e^{-z}}{z}\,dz$ is the exponential integral. Since $\{X_k\}$ are i.i.d. RVs, then the CF of the RV $Y$ is given by

$$\psi_Y(t)=\left[1+je^{-jt}E_1(-jt)\right]^K=\sum_{m=0}^K{K\choose m}j^me^{-mjt}E_1^m(-jt)$$.

Now I want to find the CDF of $Y$ from this CF. I found this theorem (Gil-Pelaez)

$$F_Y(y)=\frac{1}{2}-\frac{1}{\pi}\int_0^{\infty}\frac{\text{Im}\left\{e^{-jty}\psi_Y(t)\right\}}{t}\,dt$$

where Im(.) is the imaginary part of a complex number. However, evaluating the above integral numerically gave results that are not real and between 0 and 1, which is supposed to be. Why? What is wrong in the integral? Some said that dividing by t in the integral, while having the limits start from 0 is the problem, but this is the theorem of inverting the CF! Did I make any mistake in the derivations?

Thanks in advance

 

[mathman]

See what happens if you used this formula to get the distribution function for one X.. Also try to get the density function. Also I don''t see why the CDF is coming out not real. You should take the imaginary part, as a real function.

 

[EngWiPy]

See what happens if you used this formula to get the distribution function for one X.. Also try to get the density function. Also I don''t see why the CDF is coming out not real. You should take the imaginary part, as a real function.

When I tried to find the density of one of the Xs from the characteristic function using Mathematica, it gave me DiracDelta[x], which is wrong, because the PDF is $f_X(x)=1/(1+x)^2$. Does it mean that there is something wrong in my derivation? The derivation is shown above, and I don't see where I could have made a mistake!

 

[mathman]

I tried to check your result for char. function for one X. I found a discrepancy, specifically the si and ci terms should include the values at x=∞

 

[EngWiPy]

I tried to check your result for char. function for one X. I found a discrepancy, specifically the si and ci terms should include the values at x=∞

What do you mean? Could you elaborate, please?

 

[mathman]

I use Gradshteyn and Ryzhik as a reference. Their definite integral agrees with you, but the indefinite integral requires the extra terms when the limits are applied. I am confused myself. In any case I don't understand why you don't get the original density for one X when you take the inverse transform. If you look up the reference, the indefinite integrals are in section 2.641, wile the definite integrals are in section 3.722.

One further question - how did you get the density function of one X?

 

[EngWiPy]

... In any case I don't understand why you don't get the original density for one X when you take the inverse transform. If you look up the reference, the indefinite integrals are in section 2.641, wile the definite integrals are in section 3.722.
...

That is what I don't understand. Why I don't get the same density from the CF after inverse Fourier Transform. What reference are you talking about here?

...
One further question - how did you get the density function of one X?

$$f_X(x)=\frac{d\,F_X(x)}{d\,x}=\frac{d}{d\,x}\text{Pr}[X\leq x]=\frac{d}{d\,x}\text{Pr}[a/b\leq x]=\frac{d}{d\,x}\int_0^{\infty}\underbrace{\text{Pr}[a/\beta\leq x/b=\beta]}_{F_a(x\,\beta)}f_b(\beta)\,d\,\beta$$

If you substitute in the above equation $F_a(x)=1-e^{-x}$, and $f_b(x)=e^{-x}$, and then evaluate the integral, and take the derivative with respect to $x$ you will get $f_X(x)=1/(1+x)^2$.

 

[mathman]

Table of Integrals, Series, and Products, I.S. Gradshteyn and I.M. Ryzhik, Seventh Edition
I downloaded it a while ago.

When you did the inverse transform, did you try getting the density function? That should be easier than the CDF.
I.

 

[EngWiPy]

Table of Integrals, Series, and Products, I.S. Gradshteyn and I.M. Ryzhik, Seventh Edition
I downloaded it a while ago.

...I.

Yes, I used the definite integrals in 3.722.

...
When you did the inverse transform, did you try getting the density function? That should be easier than the CDF.
I.

Yes, when I didn't get correct answer for the CDF, I tried getting the PDF from the CF using inverse Fourier transform in Mathematica. It didn't give me the same density function!

 

[mathman]

I'll have to see more details to try to guess what happened.

 

[EngWiPy]

I think I forgot o include $t$ in the CF of $X_k$. The CF is

$$\psi_{X_k}(t)=1+jte^{-jt}E_1(-jt)$$

When I find the inverse Fourier transform of this in Mathematica as

PsiY = 1 + I t Exp[-I t] ExpIntegralE[1, -I t];
InverseFourierTransform[PsiY, t, w]

it gives me the following result

Sqrt[2 \[Pi]] DiracDelta[w]

which is not what I expect!

 

[mathman]

Have you tried using the previous equation (involving cos, sin, ci, si)? I have not used Mathematica, so I don't know what is happening there.
The Mathematica expression is for a jump distribution at 0 (after constant fixing).

 

[EngWiPy]

Have you tried using the previous equation (involving cos, sin, ci, si)? I have not used Mathematica, so I don't know what is happening there.
The Mathematica expression is for a jump distribution at 0 (after constant fixing).

I have tried this, but Mathematica ran for a while (maybe 10 minutes), and nothing was produced. My computer made a weird noise while do the computation! I would guess the CPU and memory were fully utilized for the calculation.

 

[mathman]

Based on what you are getting from Mathematica, my guess is that it is not working properly for this char. function. The Fourier transform of a delta function is a constant. Possibly Mathematica cannot handle a function of the form 1+ something. It looks to me that after the original integration by parts, you have an expression where the parts themselves are not invertable and Mathematica can't handle it.

I suggest you post your question on the following Forum. There are many more mathematicians than here who might have an answer..
https://math.stackexchange.com/

 

[EngWiPy]

I did post it on stack exchange at the same time when I posted this question here, but my question got no single response there!

 

[mathman]

I am at a loss as to what to suggest. In one of the notes, you tried to get the CDF from the char. function. Have you tried it for one X?

 

[mathman]

The expression for the Gil-Palaez form for the distribution function is real. For one X it looks like this:
$$F_X(x)=\frac{1}{2}-\frac{1}{\pi}\int_0^{\infty}(\frac{sin(tx)}{t}-si(t)sin(xt+t)-ci(t)cos(xt+t))dt$$.
See what you can do with it.

 

[mathman]

Simplification: $$\frac{1}{\pi}\int_0^\infty \frac{sin(tx)}{t}dt=\frac{1}{2}$$for all x > 0, so $$F_X(x)$$ is remaining term.

 

[EngWiPy]

The other terms seem to evaluate to simple form values depending on the arguments. For example:

$$\int_0^{\infty}\text{si}(qx)\sin(px)\,dx=\left\{\begin{array}{ll}-\frac{\pi}{2p}&p^2>q^2\\-\frac{\pi}{4p}&p^2=q^2\\0&p^2<q^2\end{array}\right.$$

Which is nice. I ill check this more thoroughly later.

 

[mathman]

Since x> 0, p(=1+x) > q(=1). Both integrals are the same, so the net result is $$\frac{-1}{1+x}$$. It looks almost like the right answer which is $$\frac{x}{1+x}$$.. This could be explained if I have the wrong sign for term leading to the integral of $$\frac{sin(t)}{t}$$. In that case, the integral becomes what it should be.

I just checked my notes and I believe I do have the wrong sign for that term. I suggest you try to get the integrand for Gil-Palaez.

 

[EngWiPy]

I will check this, but $x\geq 0$, not $x>0$. I think we need to take into account the two possibilities that $p=q$ or $p>q$, don't we?

 

[EngWiPy]

The expression for the Gil-Palaez form for the distribution function is real. For one X it looks like this:
$$F_X(x)=\frac{1}{2}-\frac{1}{\pi}\int_0^{\infty}(\frac{sin(tx)}{t}-si(t)sin(xt+t)-ci(t)cos(xt+t))dt$$.
See what you can do with it.

I checked this, and it is supposed to be like this

$$F_X(x)=\frac{1}{2}+\frac{1}{\pi}\int_0^{\infty}\left(\frac{\sin(tx)}{t}+\text{si}(t)\sin(xt+t)+\text{ci}(t)\cos(xt+t)\right)dt$$

which results in

$$F_X(x)=1-\frac{1}{x+1}$$

for $x>0$, which is the exact CDF of $X$. But if $x=0$, it yields

$$F_X(x)=1-\frac{1}{2(x+1)}$$

This is good, though. Now how can I use this insight for the summation of $X$s, instead of one $X$?

 

[mathman]

Your expression for $$F_X(x),\ x=0$$ is wrong. $$\int_0^\infty \frac{sin(tx)}{t}dt=0$$ at that point. The only insight I can give is the expansion to get the imaginary part you need will have the $$\frac{sin(tx)}{t}$$ as the first term while the rest will not be divided by t.

 

[EngWiPy]

Right. So, at $x=0$, $F_X(x)=0.5-\frac{0.5}{x+1}$, which is a scaled version of the original CDF.

 

[mathman]

$$At\ x=0, F_X(x)=0.5-0.5=0$$. I haven't looked at it in detail, but I suspect that it will also =0 for x < 0.

 

[EngWiPy]

$$At\ x=0, F_X(x)=0.5-0.5=0$$. I haven't looked at it in detail, but I suspect that it will also =0 for x < 0.

How is that $F_X(x)=0$ at $x=0$? and at $x<0$ $F_X(x)=0.5$ since

$$\int_0^{\infty}\text{si}(t)\sin(xt+t)dt=\int_0^{\infty}\text{ci}(t)\cos(xt+t)dt=0$$

for $x<0$!

 

[mathman]

$$\frac{1}{\pi}\int_0^{\infty}\frac{sin(xt)}{t}dt=-0.5,\ x\lt 0$$. You seem to have assumed = 0.

 

[EngWiPy]

OK, I see. What about the case when $x=0$. The CDF is scaled down by a factor of $0.5$, which isn't the original CDF.

 

[mathman]

Look at post #27. $$\int_0^{\infty}\frac{sin(tx)}{t}dt=0,\ for\ x=0$$

 

[EngWiPy]

Ooops! I forgot about that. My mistake. Thanks!

So, this is an indication that the derivation is correct, right? The only difference in my equations is that I use the formula $\text{ci}(x)\pm j\,\text{si}(x) = \text{Ei}(\pm j\,x)$, where $\text{Ei}(x) = -\int_{-x}^{\infty}\frac{e^{-t}}{t}\,dt = -E_1(-x)$. The question is: why the numerical integration doesn't evaluate as a CDF for the general case using MATLAB and Mathematica? Is it a precision issue in the software tools?

 

[mathman]

I don't know how MATLAB or Mathematica were set up. There may be some limitations to their applicability. Since I have never used them I can't help here. When I need complicated integrals I look up the table I previously mentioned.

 

[EngWiPy]

I don't know how MATLAB or Mathematica were set up. There may be some limitations to their applicability. Since I have never used them I can't help here. When I need complicated integrals I look up the table I previously mentioned.

But theoretically, the derivation should be correct, right? Because I just raise the CF of a single X to the power K, and then solve for the CDF as we did for a single X.

 

[mathman]

The basic theory (invert the nth power of the char. function gives the cdf of the sum of n terms) is correct. But as we have seen with the case n=1, the calculation can be very tricky. The Gil-Pelaez theorem had to be used and getting the imaginary part can be messy.