b14 change of variables plot

b14_change_of_variables.py

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

import scipy.special as special

import scipy

plt.rcParams["figure.figsize"] = (16,10)

a=861000
b=500000

v = np.linspace(a-b, a+b, 1000)
theta = np.power(4*v, 1/3)
tt = 4/3* np.power(4*v, -2/3)
t = stats.norm(loc=152.04, scale=73.5**.5).pdf(theta)
y = t * tt

# Finding optimal std for the normal curve for comparison
# fun will compute KL divergence between the pdfs: normal and `y`
fun = lambda std: stats.entropy(stats.norm(loc=v[np.argmax(t)], scale=std).pdf(v), y)
res = scipy.optimize.minimize(fun=fun, x0=[150000], method='Nelder-Mead')

ttt = stats.norm(loc=v[np.argmax(t)], scale=res.x[0]).pdf(v)

plt.plot(v, y/max(y), label='with change of variables')
plt.plot(v, t/max(t), label='without change of variables')
plt.plot(v, ttt/max(ttt), linestyle='--', label='normal curve for comparison')
plt.legend()

plt.ylabel('p(v)')
plt.xlabel('v')

plt.show()