# Tutorial: intrinsic distribution

In this tutorial you will learn:

• How to find a intrinsic distribution from data with asymmetric error bars and upper limits

• How to use UltraNest

• How to test whether there are two sub-populations

Lets say we want to find the intrinsic velocity dispersion given some noisy data points.

Our data are velocity measurements of a few globular cluster velocities in a dwarf galaxy.

:

import numpy as np
# velocity dispersions of dwarf galaxies by van Dokkum et al., Nature, 555, 629 https://arxiv.org/abs/1803.10237v1

values = np.array([15, 4, 2, 11, 1, -2, -1, -14, -39, -3])
values_lo = np.array([7, 16, 6, 3, 6, 5, 10, 6, 11, 13])
values_hi = np.array([7, 15, 8, 3, 6, 6, 10, 7, 14, 14])

n_data = len(values)


## Visualise the data

Lets plot the data first to see what is going on:

:

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure()
xlabel = 'Velocity [km/s]'
plt.xlabel(xlabel)
plt.errorbar(x=values, xerr=[values_lo, values_hi], y=range(n_data),
marker='o', ls=' ', color='orange')
plt.xlim(-50, 50); ## Data properties

This scatter plot shows:

• large, sometimes asymmetric error bars

• intrinsic scatter

### Resampling the data

We could also represent each data point by a cloud of samples. Each point represents a possible true solution of that galaxy.

:

np.random.seed(42)

:

samples = []

for i in range(n_data):
# draw normal random points
u = np.random.normal(size=400)
v = values[i] + np.where(u < 0, u * values_lo[i], u * values_hi[i])

samples.append(v)

samples = np.array(samples)

plt.figure()
# for each galaxy, plot alittle cloud with its own colors
plt.violinplot(samples.transpose(), vert=False, showextrema=False)

plt.xlabel(xlabel); ## Model

Lets fit a intrinsic, gaussian distribution.

$y \sim \mathrm{Normal}(\mu, \sigma)$

The model has two unknown parameters:

• the mean $$\mu$$

• the scatter $$\sigma$$

Lets write down prior ranges for these parameters:

:

parameters = ['mean', 'scatter']

def prior_transform(cube):
# the argument, cube, consists of values from 0 to 1
# we have to convert them to physical scales

params = cube.copy()
# let slope go from -3 to +3
lo = -100
hi = +100
params = cube * (hi - lo) + lo
# let scatter go from 1 to 1000
lo = np.log10(1)
hi = np.log10(1000)
params = 10**(cube * (hi - lo) + lo)
return params


Define the likelihood, which measures how far the data are from the model predictions. More precisely, how often the parameters would arise under the given parameters.

:

import scipy.stats

def log_likelihood(params):
# unpack the current parameters:
mean, scatter = params

# compute the probability of each sample
probs_samples = scipy.stats.norm(mean, scatter).pdf(samples)
# average over each galaxy, because we assume one of the points is the correct one (logical OR)
probs_objects = probs_samples.mean(axis=1)
assert len(probs_objects) == n_data
# multiply over the galaxies, because we assume our model holds true for all objects (logical AND)
# for numerical stability, we work in log and avoid zeros
loglike = np.log(probs_objects + 1e-100).sum()
return loglike


This likelihood evaluates the model for each object at each sample, and averages the model densities. This corresponds to saying: I don’t know which of these options is correct, but I assume one of them is, so by using a logical OR, I sum the probabilities.

The objects are then combined by multiplying the per-object probabilities. This corresponds to saying: I know the model has to describe all of these objects, so using a logical AND, I multiply the probabilities. Here we convert to log space, and sum. The small number is added to gracefully handle when a object probability is zero (distant means and small sigma).

## Solving the problem

:

import ultranest

sampler = ultranest.ReactiveNestedSampler(parameters, log_likelihood, prior_transform)


Lets first try with relatively poor sampling:

:

result = sampler.run()

[ultranest] Sampling 400 live points from prior ...
[ultranest] Widening roots to 404 live points (have 400 already) ...
[ultranest] Sampling 4 live points from prior ...

[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 5598
[ultranest]   logZ = -44.4 +- 0.07984
[ultranest] Effective samples strategy satisfied (ESS = 1720.5, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.45+-0.08 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.08, need <0.5)
[ultranest]   logZ error budget: single: 0.10 bs:0.08 tail:0.01 total:0.08 required:<0.50
[ultranest] done iterating.

:

from ultranest.plot import cornerplot
cornerplot(result) :

# ask for more samples: increasing the required effective sample size
result = sampler.run(min_ess=4000)
cornerplot(result)

[ultranest] Widening roots to 400 live points (have 404 already) ...
[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 5598
[ultranest]   logZ = -44.4 +- 0.06301
[ultranest] Effective samples strategy wants to improve: -50.96..-39.38 (ESS = 1720.5, need >4000)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.07 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.06, need <0.5)
[ultranest]   logZ error budget: single: 0.10 bs:0.06 tail:0.01 total:0.06 required:<0.50
[ultranest] Widening from 350 to 800 live points before L=-5e+01...
[ultranest] parent value is -inf, so widening roots
[ultranest] Widening roots to 800 live points (have 404 already) ...
[ultranest] Sampling 396 live points from prior ...
[ultranest] Exploring (in particular: L=-inf..-39.38) ...

[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 13079
[ultranest]   logZ = -44.32 +- 0.04902
[ultranest] Effective samples strategy wants to improve: -48.64..-39.38 (ESS = 3366.4, need >4000)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.05 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.05, need <0.5)
[ultranest]   logZ error budget: single: 0.07 bs:0.05 tail:0.00 total:0.05 required:<0.50
[ultranest] Widening from 724 to 1582 live points before L=-5e+01...
[ultranest] parent value is -inf, so widening roots
[ultranest] Widening roots to 1582 live points (have 800 already) ...
[ultranest] Sampling 782 live points from prior ...
[ultranest] Exploring (in particular: L=-inf..-39.38) ...
[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 25665
[ultranest]   logZ = -44.37 +- 0.0352
[ultranest] Effective samples strategy satisfied (ESS = 6760.5, need >4000)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.04 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.04, need <0.5)
[ultranest]   logZ error budget: single: 0.05 bs:0.04 tail:0.00 total:0.04 required:<0.50
[ultranest] done iterating. What is the 3 sigma upper limit on the scatter?

:

mean_samples, scatter_samples = result['samples'].transpose()

# get 3 sigma quantile
quantile = scipy.stats.norm().cdf(3)

# look at the value:
print('scatter is < %.4f km/s at 3 sigma (%.3f%% quantile)' % (scipy.stats.mstats.mquantiles(scatter_samples, quantile), quantile*100))

scatter is < 32.2977 km/s at 3 sigma (99.865% quantile)

:

plt.figure()

plt.subplot(2, 1, 1)
plt.errorbar(x=values, xerr=[values_lo, values_hi], y=range(n_data),
marker='o', ls=' ', color='orange')
plt.xlim(-100, 100)

plt.subplot(2, 1, 2)
plt.xlim(-100, 100)
plt.xlabel(xlabel)

from ultranest.plot import PredictionBand
x = np.linspace(-100, 100, 400)
band = PredictionBand(x)
bandc = PredictionBand(x)

for params in sampler.results['samples'][:40]:
mean, scatter = params 