Simulation of PDF of inverse random varible

Question

Let $X$ be exponential random variable with parameter $\beta$ $$ P_X(x)=\frac{1}{\beta}e^{-x/\beta}. $$

I have found a lemma saying

Let $X$ be a random variable with pdf $p_X (x)$ for all $x\geq 0$ and $p_X (x) = 0$ for $x < 0$. Then, the pdf of $Y = 1/ X$ is $$ P_Y(y)=\frac{1}{y^2}P_X(\frac{1}{y}). $$

My question is, How we use MATLAB with Monte Carlo simulation to generate or simulate $X$ and compared with theoretical PDF $P_X(x)$ and also simulate $P_Y(y)$ and compared with theoretical expression.

Answer 1

The classical way to simulate an $Exp[\beta]$ random variable is by considering random variable

$$X'=-\beta \log(U)$$

where $U$ is uniformly distributed on $[0,1] \ $ ($\log$ is the natural logarithm). It is easy to prove that $X'$ has a $Exp[\beta]$ distribution.

Thus the simulation program in Matlab is :

 clear all;close all;
 n=100000;
 beta=10;
 E=-beta*log(rand(1,n));
 hist(E,30)

---

Second part : What can be said about $Z$ defined by

$$\dfrac{1}{Z}=\dfrac{1}{X}+\dfrac{1}{Y}$$

with $X,Y$ following a $Exp[\beta]$ distribution ?

The problem with the inverse of an $Exp[\beta]$ distribution (like $\dfrac{1}{X}$ or $\dfrac{1}{Y}$) is that it has an infinite mean, thus a non-defined variance (https://stats.stackexchange.com/questions/229543/mean-of-inverse-exponential-distribution)

Thus, instead of simulating $1/X$; $1/Y$, adding them and taking the inverse, I have directly done this simulation, for the particular case $\beta=1$ :

 clear all;close all;
 n=1000000;
 U=-log(rand(1,n));V=-log(rand(1,n));
 Z=(U.*V)./(U+V);
 hist(Z,200)
 [mean(Z),var(Z)]

The results I obtain tend to give $E(Z)=\dfrac{1}{3}$ and $var(X)=\dfrac{4}{45}$.

Here is the histogram generated by this program. One notices that after a certain small peak, there is rapid decline.

I have no idea about the way to get these results throretically.