In [92]:

    

import sympy as sy
from IPython.display import display
sy.init_printing(use_latex='mathjax')


    

In [93]:

    

s, lambda_ = sy.symbols('s lambda', real=True, positive=True)
beta = sy.symbols('beta', real=True)


    

In [94]:

    

gamma = lambda_**2/2 * sy.exp(-lambda_**2 / 2 * s)
display(gamma)


    

    






$$\frac{\lambda^{2}}{2 e^{\frac{\lambda^{2} s}{2}}}$$

In [95]:

    

normal = 1/s*sy.exp(-beta**2/2/s**2)
display(normal)


    

    






$$\frac{1}{s} e^{- \frac{\beta^{2}}{2 s^{2}}}$$

In [96]:

    

sy.simplify(normal * gamma)


    

    
Out[96]:





$$\frac{\lambda^{2}}{2 s e^{\frac{\beta^{2}}{2 s^{2}} + \frac{\lambda^{2} s}{2}}}$$

In [97]:

    

sy.integrate(normal * gamma, s)


    

    
Out[97]:





$$\int \frac{\lambda^{2} e^{- \frac{\beta^{2}}{2 s^{2}}}}{2 s e^{\frac{\lambda^{2} s}{2}}}\, ds$$

In [98]:

    

sy.expand_log(2 * sy.log(normal * gamma))


    

    
Out[98]:





$$- \frac{\beta^{2}}{s^{2}} - \lambda^{2} s + 4 \log{\left (\lambda \right )} - 2 \log{\left (s \right )} - 2 \log{\left (2 \right )}$$

In [99]:

    

%matplotlib inline
import scipy.stats as st
import matplotlib.pyplot as plt
import numpy as np
N = 1000000
r1 = st.norm(0, 1).rvs(N)
r2 = st.gamma(.5, .5).rvs(N)
r = st.expon(1).rvs(N)


    

In [100]:

    

_ = plt.hist(r1)


    

In [101]:

    

_ = plt.hist(r2)


    

In [102]:

    

_ = plt.hist(r)


    

In [106]:

    

_ = plt.hist(r1 * np.sqrt(r/r2), bins=1000)
plt.xlim(-1.5, 1.5)


    

    
Out[106]:





$$\begin{pmatrix}-1.5, & 1.5\end{pmatrix}$$






    

In [ ]: