```python
%pylab
```

    Using matplotlib backend: Qt5Agg
    Populating the interactive namespace from numpy and matplotlib



```python
from random import uniform
```

Given $X \sim \Gamma(\alpha,1)$, with density
\begin{equation}
f_X(x) \propto f^*_X(x) = x^{\alpha-1} e^{-x}, \quad x > 0
\end{equation}
Using the envelope $Y \sim \text{Exp}(1/\alpha)$ with density
\begin{equation}
g_Y(x) = \frac{1}{\alpha} e^{-x/\alpha}, x > 0
\end{equation}
Also compute the probability of rejection.


```python
f = lambda alpha,x: (x**(alpha-1))*(np.exp(-x))
g = lambda alpha,x: (1/alpha)*(np.exp(-x/alpha))
```

1. Compute $\displaystyle{ M = \sup_{x\in(0,1)} \frac{f^*_X(x)}{g_Y(x)}}$: <br>
Note that 
\begin{equation}
h(x) = \frac{f^*_X(x)}{g_Y(x)} = \frac{1}{\alpha}x^{\alpha-1}\exp\left(-x\left(1-\frac{1}{\alpha}\right)\right)
\end{equation}
To find supremum we may consider the quantity
\begin{equation}
y = \ln h(x) = - \ln \alpha + (\alpha - 1) \ln x - x \left( 1-\frac{1}{\alpha} \right)
\end{equation}
Then 
\begin{equation}
\frac{dy}{dx} = \frac{\alpha - 1}{x} - \frac{\alpha-1}{\alpha}
\end{equation}
This is equal to zero if $\displaystyle{x = \alpha}$. In fact, 
    - if $x < \alpha$ then $dy/dx < 0$
    - if $x > \alpha$ then $dy/dx > 0$
Therefore $y$ (and hence $h$) attains maximum at $x = \alpha$, therefore $\displaystyle{M = \alpha^{\alpha-2}e^{-(\alpha-1)}}$

2. We may simulate Y by using inversion - 
Note that the CDF of $Y$ is $G_Y(x) = 1 - e^{-x/\alpha}$, and $G_Y^{-1}(x) = -\alpha \ln(1-x)$


```python
def my_rExp(alpha,size):
    unisamp = [random.uniform(0,1) for i in range(size)]
    return [-alpha*np.log(1-x) for x in unisamp]
```

3. We may then use rejection algorithm to simulate $X$.
    - Simulate $U \sim U(0,1)$ and $X \sim \text{Exp}(1/\alpha)$.
    - If $u \leq M^{-1}f^*_X(x)/g(x)$ then accept $u$, otherwise reject.


```python
M = lambda alpha: (alpha**(alpha-2))*np.exp(-(alpha-1))
h = lambda alpha,x: f(alpha,x)/(M(alpha)*g(alpha,x))
```


```python
def my_rGamma1(alpha,size):
    count = 0
    samp = []
    while count < size:
        u = [random.uniform(0,1) for i in range(size)]
        exp_samp = my_rExp(alpha,size)
        ubound = [h(alpha,x) for x in exp_samp]
        idx = [i for i in range(len(u)) if u[i] <= ubound[i]]
        new_samp = [u[i] for i in idx]
        samp += new_samp
        count = len(samp)
    return [samp[i] for i in range(size)]
```


```python
my_rGamma1(2,10)
```




    [0.5550633830391315,
     0.9927212674394112,
     0.6983735483799777,
     0.49466601339924643,
     0.8106875849861941,
     0.9050509249637454,
     0.3102623813753336,
     0.98171935773631,
     0.19668951612576313,
     0.36127920327731633]



Notice that if $\tilde{X} = X/\beta$ then $\tilde{X} \sim \Gamma(\alpha,\beta)$. Using scaling we may simulate $\tilde{X}$.


```python
def my_rGamma(alpha,beta,size):
    return [pt/beta for pt in my_rGamma1(alpha,size)]
```


```python
my_rGamma(3,2,30)
```




    [0.38359214592412577,
     0.4966073373612345,
     0.2568435931525481,
     0.3975477716675325,
     0.31461628807870534,
     0.43438550365391143,
     0.05458809067372872,
     0.33602126523235437,
     0.22906344692827613,
     0.3177055644399481,
     0.21831892032313138,
     0.3036546829741233,
     0.16840905801462613,
     0.15654361943179157,
     0.2767799227870573,
     0.1533971950751764,
     0.04082917886722137,
     0.07897936636851083,
     0.4179233869683552,
     0.2489387208380447,
     0.31459121925637584,
     0.002321617317671887,
     0.028271625135758538,
     0.20481184134471248,
     0.32464579215021194,
     0.1317549529978439,
     0.23665240812238492,
     0.44853126271481225,
     0.4848416365454395,
     0.39500027839659013]



02f6c0c9bcf0b18835ec94c9e94b513019b3c642
ipynb
Jupyter Notebook
M3S9/Rejection Sampling Gamma(a,b) Example.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
2021-05-16T04:08:34.000Z
2021-11-27T12:56:10.000Z
M3S9/Rejection Sampling Gamma(a,b) Example.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
M3S9/Rejection Sampling Gamma(a,b) Example.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
2020-03-31T00:23:13.000Z
2020-05-13T15:01:46.000Z
Qwen/Qwen-72B
1. YES
2. YES