<jupyter_start><jupyter_text>Example 0: Fitting a line to data --- 1D Gauss likelihood, flat priorAs the first example, lets take the simplest possible scenario of fitting a line to data. $${\rm data} = \{(x_i, y_i)\}_{i = 1}^{N_{\rm data}}$$where $y$ has a linear relation with $x$ with some constant Gaussian noise$$y \sim \mathcal{N}(-x + 0.5, \sigma)$$ We'll fix $\sigma = 0.1$. Meanwhile, our model will be equally simplistic: $$m(x, \theta) = -x + \theta$$So we're just inferring the $y$-intercept of a line. Then the likelihood is a 1D Gaussian: $$p(d|\theta) = \mathcal{L}(\theta) = \prod\limits_{i=1}^{N_{\rm data}}\frac{1}{\sqrt{2\pi\sigma_i^2}} exp\left(-\frac{(d_i - m(\theta)_i)^2}{2\sigma_i^2} \right)$$and we'll use a simple flat prior $$p(\theta) = 1~{\rm if}~0 < \theta < 1; 0~{\rm else}$$With this flat prior, things become pretty trivial. The evidence, which is the integral$$p({\rm data}) = \int_{\Omega_\theta} p({\rm data}|\theta) p(\theta)~{\rm d}\theta$$($\Omega_\theta$ is the parameter space defined by the flat prior), becomes$$p({\rm data}) = \int_{\Omega_\theta} p({\rm data}|\theta)~{\rm d}\theta$$Things should be equally trivial for the log evidence estimator, since the KL divergence in this case simply becomes$$D_{\rm KL} (p(\theta|{\rm data} || p(\theta) = \int_{\Omega_\theta} p(\theta |{\rm data}) \log \frac{p(\theta |{\rm data})}{p(\theta)} {\rm d}\theta = \int_{\Omega_\theta} p(\theta |{\rm data}) \log p(\theta |{\rm data}) {\rm d}\theta$$<jupyter_code>theta_true = 0.5

def model(theta, x): 
    '''forward model
    '''
    return -x + theta

N_data = 25

np.random.seed(0)
sig_noise = 0.2
x_data = 10.*np.random.uniform(size=N_data)-5.
y_data =  model(theta_true, x_data) + sig_noise * np.random.randn(N_data)
y_err = np.repeat(sig_noise, N_data)

fig = plt.figure(figsize=(10,5))
sub = fig.add_subplot(111)
sub.plot(x_data, model(theta_true, x_data), c='C0')
sub.errorbar(x_data, y_data, yerr=y_err, fmt='.k')
sub.set_xlabel('x', fontsize=25)
sub.set_xlim(-5, 5)
sub.set_ylabel('y', fontsize=25)
prior_lim = [[0., 1.]]
def lnlikelihood(theta): 
    ''' gaussian likelihood
    '''
    y_model = model(theta, x_data)
    chi2 = np.sum((y_data - y_model)**2/y_err**2)
    
    lnlike = -0.5 *chi2 + np.log(2.*np.pi*y_err**2).sum()
    return lnlike

def likelihood(theta): 
    return np.exp(lnlikelihood(theta))

def lnprior(theta): 
    '''flat prior
    '''    
    if prior_lim[0][0] < theta < prior_lim[0][1]:
        _prior = 1./(prior_lim[0][1] - prior_lim[0][0])
        return np.log(_prior)
    return -np.inf

def lnposterior(theta): 
    return lnlikelihood(theta) + lnprior(theta)
N_prior = 10000
prior_sample = np.array([np.random.uniform(lim[0], lim[1], size=N_prior) for lim in prior_lim]).T
fig = corner.corner(prior_sample, labels=['$b$'])
from scipy.optimize import minimize

nll = lambda *args: -lnlikelihood(*args)
initial = np.array(theta_true) + 0.1 * np.random.randn(1)

soln = minimize(nll, initial)
theta_ml = soln.x

pos = soln.x + 1e-4 * np.random.randn(32, 1)
nwalkers, ndim = pos.shape

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnposterior)
sampler.run_mcmc(pos, 5000, progress=True)
chain = sampler.get_chain(discard=100, thin=10, flat=True)
fig = corner.corner(chain, labels=['$b$'], truths=[theta_true])
lnlike_chain = np.array([lnlikelihood(tt) for tt in chain])<jupyter_output><empty_output><jupyter_text>Lets first calculate the evidence by explicitly evaluating the integral$$p({\rm data}) = \int_{\Omega_\theta} p({\rm data}|\theta) p(\theta)~{\rm d}\theta$$which in our case is $$p({\rm data}) = \int_{\Omega_\theta} p({\rm data}|\theta)~{\rm d}\theta$$<jupyter_code>evidence, _ = integrate.quad(likelihood, prior_lim[0][0], prior_lim[0][1])
evidence /= (prior_lim[0][1] - prior_lim[0][0])
print('evidence from integrating the likelihood')
print(evidence)
uniq_chain, iuniq = np.unique(chain, return_index=True)

log_evid = Nevin.log_evidence(np.atleast_2d(uniq_chain).T, lnlike_chain[iuniq], prior_sample, silent=False)
print('evidence from estimator')
print(np.exp(log_evid))<jupyter_output><ln L> = -45.704759
  epsilons  [0.01278185 0.00031637 0.00031637 ... 0.00062836 0.00288296 0.00648057]
  l_i  [1 1 1 ... 1 1 1]
  k_i  [255  14  12 ...  16  57 113]
  first term = 0.115073
  digamma term = 2.132672
D_KL = 1.798071
evidence from estimator
2.3430585561754e-21<jupyter_text>**estimator is consistent**.  Example 2: Fitting a line to data --- 2D Gauss likelihood, flat priorsame as above but now slope is a free parameter<jupyter_code>theta_true = [-1, 0.5]

def model(theta, x): 
    '''forward model
    '''
    return theta[0] * x + theta[1]

N_data = 25

np.random.seed(0)
sig_noise = 0.2
x_data = 10.*np.random.uniform(size=N_data)-5.
y_data =  model(theta_true, x_data) + sig_noise * np.random.randn(N_data)
y_err = np.repeat(sig_noise, N_data)

fig = plt.figure(figsize=(10,5))
sub = fig.add_subplot(111)
sub.plot(x_data, model(theta_true, x_data), c='C0')
sub.errorbar(x_data, y_data, yerr=y_err, fmt='.k')
sub.set_xlabel('x', fontsize=25)
sub.set_xlim(-5, 5)
sub.set_ylabel('y', fontsize=25)
prior_lim = [[-5, 0], [0., 1.]]
def lnlikelihood(theta): 
    ''' gaussian likelihood
    '''
    N = theta.shape[0]
    y_model = model(theta, x_data)
    chi2 = np.sum((y_data - y_model)**2/y_err**2)
    lnlike = -0.5 * chi2 + np.log(2.*np.pi*y_err**2).sum()
    return lnlike

def likelihood(theta): 
    return np.exp(lnlikelihood(theta))

def lnprior(theta): 
    '''flat prior
    '''    
    if prior_lim[0][0] < theta[0] < prior_lim[0][1] and prior_lim[1][0] < theta[1] < prior_lim[1][1]:
        _prior = 1./(prior_lim[0][1] - prior_lim[0][0])/(prior_lim[1][1] - prior_lim[1][0])
        return np.log(_prior)
    return -np.inf

def lnposterior(theta): 
    return lnlikelihood(theta) + lnprior(theta)
N_prior = 10000
prior_sample = np.array([np.random.uniform(lim[0], lim[1], size=N_prior) for lim in prior_lim]).T
fig = corner.corner(prior_sample, labels=['$m$', '$b$'])
nll = lambda *args: -lnlikelihood(*args)
initial = np.array(theta_true) + 0.1 * np.random.randn(2)

soln = minimize(nll, initial)
theta_ml = soln.x

pos = soln.x + 1e-4 * np.random.randn(32, 2)
nwalkers, ndim = pos.shape

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnposterior)
sampler.run_mcmc(pos, 5000, progress=True)
chain = sampler.get_chain(discard=100, thin=10, flat=True)
fig = corner.corner(chain, labels=['$m$', '$b$'], truths=theta_true)
_likelihood = lambda y, x: likelihood(np.array([y,x]))
evidence, _ = integrate.dblquad(_likelihood, prior_lim[1][0], prior_lim[1][1], lambda x: prior_lim[0][0], lambda x: prior_lim[0][1])
evidence /= (prior_lim[0][1] - prior_lim[0][0])
print('evidence from integrating the likelihood')
print(evidence)
lnlike_chain = np.array([lnlikelihood(tt) for tt in chain])
log_evid = Nevin.log_evidence(chain, lnlike_chain, prior_sample, silent=False)
print('evidence from estimator')
print(np.exp(log_evid))<jupyter_output><ln L> = -45.721869
  epsilons  [0.01684253 0.00721258 0.01220292 ... 0.00773628 0.00699415 0.01710189]
  l_i  [1978   99  489 ...  316  376 2208]
  k_i  [1 1 1 ... 1 1 1]
  first term = 0.008425
  digamma term = 6.896488
D_KL = 6.455176
evidence from estimator
2.1867408814599582e-23
{'fraction_non_alphanumeric': 0.17724319093045762, 'fraction_numerical': 0.060279275542651736, 'mean_word_length': 4.4803030303030305, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 13, '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': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '129896', 'license': 'MIT', 'n_tokens_mistral': 3143, 'n_tokens_neox': 2797, 'n_words': 648, 'path': 'nb/eg_easy.ipynb', 'repo_name': 'changhoonhahn/nevin', 'size': '7233'}