The PDF of random variable $Z$ is given by

$$f_{Z}(z)=\hspace{-2mm}\int \limits_0^{\infty}\frac{z^{2m-1}\exp\left(-\frac{mz^2}{\sigma}\right)}{\Gamma(m)\sigma^m}\frac{2m^m}{\sqrt{2\pi}\zeta\sigma}\exp\left(-\frac{(\ln \sigma-\mu)^2}{2\zeta^2}\right)d\sigma$$.

Let $m=2$, $\mu=0.6931$ and $\zeta=0.9240$. How can I plot this PDF?

I tried this:

μ = 0.6931;
ζ = 0.9240;
m = 2;

f[z_] :=
    (z^(2*m - 1)*Exp[-m*z^2/σ]*(2*m^m)*
      Exp[-(Log[σ] - μ)^2/(2*ζ^2)]/(Gamma[m]*σ^m*Sqrt[2*π]*ζ*σ)), 
    {σ, 0, Infinity}]

Dist = ProbabilityDistribution[f[z], {z, 0, Infinity}];

pdfH[z] := PDF[Dist, z];

Plot[pdfH, {z, 0, 5}];

It does not work.

Furthermore, if $X$ is a random variable which is given by $X = Z^2$, then how to draw the PDF of $X$?

  • $\begingroup$ How would you plot anything? Your syntax of Plot is wrong, you suppress the output of the command (the actual plot), and you are generally confusing = and := definitions. Please start by reading the manual reference.wolfram.com/language/ref/Plot.html?q=Plot $\endgroup$ – Felix Jan 26 '17 at 1:34
  • 3
    $\begingroup$ I don't think this is a well posed PDF. Doing the integration of the PDF numerically, f[z_?NumericQ] := NIntegrate[(z^(2*m - 1)*Exp[-m*z^2/\[Sigma]]*2*m^m* Exp[-(Log[\[Sigma]] - \[Mu])^2/(2*\[Sigma]^2)]/(Gamma[m]*\[Sigma]^ m*Sqrt[2*\[Pi]]*\[Zeta]*\[Sigma])), {\[Sigma], 0, Infinity}], allows for direct evaluation for particular z. However, for the PDF to be valid, the integral of the PDF over the domain of z must be equal to one. Here assuming $z \in \left[0,\infty\right]$, that is not the case. For the integral of the PDF to be 1, the upper bound of z is 2.5178. $\endgroup$ – Marchi Jan 26 '17 at 2:24
  • $\begingroup$ @Marchi, how do you find the upper bound of z? Would you please show me the way. If we fix the upper limit of z to the above value, is it possible to draw the PDF then? $\endgroup$ – Dimitrios Jan 26 '17 at 2:38
  • $\begingroup$ Your Mathematica code expressing the integral does not correspond to your MathJax expression giving the integral. $\endgroup$ – m_goldberg Jan 26 '17 at 2:59
  • $\begingroup$ @m_goldberg, I think it is correct now. $\endgroup$ – Dimitrios Jan 26 '17 at 3:06

Irrespective of the discussion if this is a well defined PDF, the issue why you don't get a plot is because the Integrate in the definition of f can not be calculated analytically. You need a numerical integral, as already suggested by Marchi:

μ = 0.6931;
ζ = 0.9240;
m = 2;

f[z_?NumericQ] := 
 NIntegrate[(z^(2*m - 1)*Exp[-m*z^2/σ]*(2*m^m)*
    Exp[-(Log[σ] - μ)^2/(2*ζ^2)]/(Gamma[m]*σ^
        m*Sqrt[2*π]*ζ*σ)), {σ, 0, Infinity}]

Dist = ProbabilityDistribution[f[z], {z, 0, Infinity}];

pdfH[z_] := PDF[Dist, z];
Plot[pdfH[z], {z, 0, 5}]

Note the z in Dist is not the same as the z in Plot. Therefore, using pattern matching and SetDelayed are important to not confuse the two.

enter image description here

  • $\begingroup$ Thank you very much. However, I get a blank plot with the above Mathematica code.. $\endgroup$ – Dimitrios Jan 26 '17 at 4:42
  • $\begingroup$ Did you use a fresh kernel? Or at least after a ClearAll? $\endgroup$ – Felix Jan 26 '17 at 4:59
  • $\begingroup$ Thank you very much. I am getting it now. I really appreciate it very much. $\endgroup$ – Dimitrios Jan 26 '17 at 5:05
  • $\begingroup$ Shouldn't this work? NExpectation[x, x \[Distributed] Dist] ? Given that this does : NIntegrate[x pdfH[x] , {x, 0, Infinity}] $\endgroup$ – chris Dec 5 '18 at 8:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.