<jupyter_start><jupyter_text>Ellipsoidal nested samplingThis example demonstrates how to use ellipsoidal nested rejection sampling [1] to sample from the posterior distribution for a logistic model fitted to model-simulated data.[1] "A nested sampling algorithm for cosmological model selection", <NAME>, <NAME> and <NAME>, [arXiv:astro-ph/0508461v2](https://arxiv.org/abs/astro-ph/0508461). First create fake data.<jupyter_code>import pints
import pints.toy as toy
import numpy as np
import matplotlib.pyplot as plt

# Load a forward model
model = toy.LogisticModel()

# Create some toy data
r = 0.015
k = 500
real_parameters = [r, k]
times = np.linspace(0, 1000, 100)
signal_values = model.simulate(real_parameters, times)

# Add independent Gaussian noise
sigma = 10
observed_values = signal_values + pints.noise.independent(sigma, signal_values.shape)

# Plot
plt.plot(times,signal_values,label = 'signal')
plt.plot(times,observed_values,label = 'observed')
plt.xlabel('Time')
plt.ylabel('Values')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Create the nested sampler that will be used to sample from the posterior.<jupyter_code># Create an object with links to the model and time series
problem = pints.SingleOutputProblem(model, times, observed_values)

# Create a log-likelihood function (adds an extra parameter!)
log_likelihood = pints.GaussianLogLikelihood(problem)

# Create a uniform prior over both the parameters and the new noise variable
log_prior = pints.UniformLogPrior(
    [0.01, 400, sigma * 0.5],
    [0.02, 600, sigma * 1.5])

# Create a nested ellipsoidal rejectection sampler
sampler = pints.NestedController(log_likelihood, log_prior, method=pints.NestedEllipsoidSampler)

# Set number of iterations
sampler.set_iterations(8000)

# Set the number of posterior samples to generate
sampler.set_n_posterior_samples(1600)

# Do proposals in parallel
sampler.set_parallel(True)

# Use dynamic enlargement factor
sampler._sampler.set_dynamic_enlargement_factor(1)<jupyter_output><empty_output><jupyter_text>Run the sampler!<jupyter_code>samples = sampler.run()
print('Done!')<jupyter_output>Running Nested ellipsoidal sampler
Number of active points: 400
Total number of iterations: 8000
Total number of posterior samples: 1600
Iter. Eval. Time m:s Delta_log(z) Acceptance rate
0     1       0:05.6 -inf          1             
0     2       0:05.6 -inf          1             
0     21      0:05.6 -inf          1             
0     41      0:05.7 -inf          1             
0     61      0:05.7 -inf          1             
0     81      0:05.8 -inf          1             
0     101     0:05.8 -inf          1             
0     121     0:05.8 -inf          1             
0     141     0:05.9 -inf          1             
0     161     0:05.9 -inf          1             
0     181     0:06.0 -inf          1             
0     201     0:06.0 -inf          1             
0     221     0:06.1 -inf          1             
0     241     0:06.1 -inf          1             
0     261     0:06.2 -inf          1             
0     281     0:06.2 -inf          1             
0     301    [...]<jupyter_text>Plot posterior samples versus true parameter values (dashed lines)<jupyter_code># Plot output
import pints.plot

pints.plot.histogram([samples], ref_parameters=[r, k, sigma])

plt.show()<jupyter_output><empty_output><jupyter_text>Plot posterior predictive simulations versus the observed data<jupyter_code>pints.plot.series(samples[:100], problem)
plt.show()<jupyter_output><empty_output><jupyter_text>Marginal likelihood estimate<jupyter_code>print('marginal log-likelihood = ' + str(sampler.marginal_log_likelihood())
      + ' ± ' + str(sampler.marginal_log_likelihood_standard_deviation()))<jupyter_output>marginal log-likelihood = -375.348895962891 ± 0.07908966826933428<jupyter_text>Effective sample size<jupyter_code>print('effective sample size = ' + str(sampler.effective_sample_size()))<jupyter_output>effective sample size = 1597.4248390822784
StarcoderdataJupyterStructuredCleanDedup