<jupyter_start><jupyter_text>Transformation of likelihood with change of random variableIn this notebook we will be looking at how the likelihood and likelihood ratio transform under a change of variables for the random variable. In particular we will consider $X\sim P_X(X|\sigma) := N(\mu=1,\sigma)$ and $Q = 1/X$. We will also use automatic differentiation to calculate $P_Q(Q|\sigma)$. by <NAME>, September 16, 2020<jupyter_code>import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams.update({'font.size': 15})
#rcParams.update({'font.size': 18,'text.usetex': True})
rcParams['figure.figsize'] = 12,6
true_sigma=0.3
true_mean = 1.5
N = 10000
def my_family_of_pdfs_for_X(x, sigma):
    #return 1./np.sqrt(2.*np.pi)/sigma * np.exp( -(x-mu)**2/2/sigma**2 )
    return norm.pdf(x,loc=true_mean, scale=sigma)
# Generate a lot of synthetic data for visualizing distribution 
x_for_hist = np.random.normal(true_mean, true_sigma, N)
# Generate a small synthetic data set
x_for_data = norm.rvs(loc=true_mean, scale=true_sigma, size=15)
x_for_data
x_for_plot = np.linspace(0.1,3,100)
likelihood_i = my_family_of_pdfs_for_X(x_for_data, true_sigma)
_ = plt.hist(x_for_hist, bins=50, density=True, alpha=.3)
plt.plot(x_for_plot, my_family_of_pdfs_for_X(x_for_plot, true_sigma), c='r', label=r'$p(X|\sigma)$',lw=2)
plt.scatter(x_for_data, likelihood_i, c='r', label=r'$L_i(\sigma) = p(x_i \mid \sigma=\sigma_{true})$')
plt.plot(x_for_data, [0.01]*len(x_for_data), '|', color='k', label=r'$x_i$')
plt.xlabel(r'$x$')
plt.ylabel(r'$p(x \mid \sigma=\sigma_{true}$)')
plt.legend()
def q(x):
    return 1/x
q_for_hist = q(x_for_hist)
q_for_plot = q(x_for_plot)
q_for_data = q(x_for_data)
from jax import grad, vmap
import jax.numpy as np
#define the gradient with grad(q)
dq = vmap(grad(q))  #dq is a new python function
print(dq(np.array([.5]))) # should be -4
def my_family_of_pdfs_for_Q(q, sigma):
    #use change of variables formua
    x = 1/q
    base_pdf = my_family_of_pdfs_for_X(x, sigma)
    jacobian = np.abs(dq(x))
    return base_pdf/jacobian
import seaborn as sns
a_plot = sns.jointplot(x=x_for_hist, y=q_for_hist, kind="hist", fill=True)
a_plot.ax_marg_y.set_ylim(0.25, 1.25)
a_plot.set_axis_labels('X', 'Q', fontsize=16)
_ = plt.hist(q_for_hist, bins=np.linspace(0,2,50), density=True, alpha=.3)
plt.plot(q_for_plot, my_family_of_pdfs_for_Q(q_for_plot, true_sigma), c='r', lw=2, label=r'$p(Q|\sigma)$')
plt.scatter(q_for_data, my_family_of_pdfs_for_Q(q_for_data, true_sigma), c='r', lw=2, label=r'$L^\prime_i(\sigma) = p(q_i \mid \sigma=\sigma_{true})$')
plt.plot(q_for_data, [0.01]*len(q_for_data), '|', color='k', label=r'$q_i$')
plt.xlim((0,2))
plt.xlabel('q')
plt.ylabel('p(q)')
plt.legend()
def my_X_likelihood(data,sigma):
    L = 1.
    for x in data:
        L *= my_family_of_pdfs_for_X(x,sigma)
    return L

'''
import numpy as np
def my_Q_likelihood(data,sigma):
    """This is a plane native python way to return array of likelihood values for each sigma"""
    ret = 0.*sigma
    for i, s in enumerate(sigma):
        ret[i] = np.prod(my_family_of_pdfs_for_Q(data,s))
    return ret
'''


#import jax.numpy as np
def my_Q_likelihood(data,sigma):
    """this is a faster way of acheiving the above with broadcasting"""
    data_ = np.tile(data, sigma.size)
    sigma_ = np.repeat(sigma, data.size)
    return np.prod(my_family_of_pdfs_for_Q(data_,sigma_).reshape(sigma.size,data.size),axis=1)
sigma_for_plot = np.linspace(.1,1)
x_likelihood = my_X_likelihood(x_for_data, sigma_for_plot)
q_likelihood = my_Q_likelihood(q_for_data, sigma_for_plot)
likelihood_offset = -2.*np.log(x_likelihood)[0]+2.*np.log(q_likelihood)[0]
print('observed offset is ', likelihood_offset)
def predict_offset_from_jacobian(data):
    return 2*np.sum(np.log(1./np.abs(dq(data))))

print("predicted offset is ", predict_offset_from_jacobian(x_for_data))
fig, ax = plt.subplots(constrained_layout=True)

ax.plot(sigma_for_plot, -2.*np.log(x_likelihood), label=r'$-2 \log L(\sigma))$', lw=2)
ax.plot(sigma_for_plot, -2.*np.log(q_likelihood), label=r'$-2 \log L^\prime(\sigma)$', ls='dashed',lw=2)
ax.plot(sigma_for_plot, np.ones(sigma_for_plot.size)*likelihood_offset, label=r'offset', ls='dashed',lw=2)
ax.set_ylabel(r'$-2 \log L$')
ax.set_xlabel(r'$\sigma$')
ax.legend()
fig, ax = plt.subplots(constrained_layout=True)

ax.plot(sigma_for_plot, -2.*np.log(x_likelihood), label=r'$-2 \log L(\sigma))$', lw=2)
ax.plot(sigma_for_plot, -2.*np.log(q_likelihood)+likelihood_offset, label=r'$-2 \log L^\prime(\sigma)$+offset', ls='dashed',lw=2)
ax.set_ylabel(r'$-2 \log L$')
ax.set_xlabel(r'$\sigma$')
ax.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1652452025586354, 'fraction_numerical': 0.020042643923240937, 'mean_word_length': 5.878299120234604, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '195435', 'license': 'MIT', 'n_tokens_mistral': 2009, 'n_tokens_neox': 1921, 'n_words': 377, 'path': 'book/distributions/likelihood-change-obs.ipynb', 'repo_name': 'wjnbreu/stats-ds-book', 'size': '4696'}