mc-variance.py

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


def sf_rad(mm, aa, bb):
    rs = aa + bb * (np.log10(mm) - 10.0)
    return np.power(10.0, rs)


# Number of Monte-Carlo iterations
NUM = int(1e3)

# These are fit parameters, with the average and +-sigma values
aa = np.array([0.171, 0.008])
bb = np.array([0.226, 0.022])

# x-values (independent variable)
mstar = np.logspace(8, 12, 40)

# Draw `NUM` variables from normal distributions for each parameter
ar = np.random.normal(aa[0], aa[1], size=NUM)
br = np.random.normal(bb[0], bb[1], size=NUM)

# Construct `NUM` curves (axis=1) for each x-value (axis=0)
# Shape (X, NUM)
rs = sf_rad(mstar[:, np.newaxis], ar[np.newaxis, :], br[np.newaxis, :])

# Median/average relationship for plotting
rs_med = sf_rad(mstar, aa[0], bb[0])

# Create a nice figure
fig, ax = plt.subplots(figsize=[10, 6])
ax.grid(True)
ax.set(xscale='log', yscale='log')
ax.tick_params(axis='x', which='both', colors='k', labelsize=14)
ax.tick_params(axis='y', which='both', colors='b', labelsize=14)
ax.set_ylabel('Values ($y$)', fontsize=20, color='b')

# Plot median relationship
ax.plot(mstar, rs_med, 'k--', lw=2.0)

# Plot every `SKIP`th MC curve (might be a little slow; make sure that the axis=0 shapes agree)
SKIP = 30
ax.plot(mstar, rs[:, ::SKIP], 'b-', alpha=0.1)

# Normal-distribution CDF as a function of std/sigma values
# i.e. this is the fraction of a normal distribution below each sigma value
fracs = sp.stats.norm.cdf([-2, -1, 0, 1, 2])

# Calculate the percentiles along the Monte-Carlo'd axis (axis=1)
percs = np.percentile(rs, 100*fracs, axis=1)
# Plot the region between -1,+1 sigma
ax.fill_between(mstar, percs[1, :], percs[3, :], color='r', alpha=0.5)

# Calculate standard-deviation (symmetric)
std = np.std(rs, axis=1)

# Plot std on a twin-y axis (normalize to y-values --- i.e. fractional-std)
tw = ax.twinx()
tw.plot(mstar, std/rs_med, 'r--')
tw.tick_params(axis='y', which='both', colors='r', labelsize=14)
tw.set_ylabel('Std / Values ($\sigma_y/y$)', fontsize=20, color='r')

plt.show()