import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import statsmodels.api as sm
# R-like interface, alternatively you can import statsmodels as import statsmodels.api as sm
import statsmodels.formula.api as smf
import theano
from scipy import stats
from scipy.special import logsumexp
%config InlineBackend.figure_format = 'retina'
az.style.use('arviz-darkgrid')
data = {'species' : ['afarensis', 'africanus', 'habilis', 'boisei', 'rudolfensis', 'ergaster', 'sapiens'],
'brain' : [438, 452, 612, 521, 752, 871, 1350],
'mass' : [37., 35.5, 34.5, 41.5, 55.5, 61.0, 53.5]}
d = pd.DataFrame(data)
d
species | brain | mass | |
---|---|---|---|
0 | afarensis | 438 | 37.0 |
1 | africanus | 452 | 35.5 |
2 | habilis | 612 | 34.5 |
3 | boisei | 521 | 41.5 |
4 | rudolfensis | 752 | 55.5 |
5 | ergaster | 871 | 61.0 |
6 | sapiens | 1350 | 53.5 |
m_6_1 = smf.ols('brain ~ mass', data=d).fit()
1 - m_6_1.resid.var()/d.brain.var()
# m_6_1.summary() check the value for R-squared
0.490158047949084
m_6_2 = smf.ols('brain ~ mass + I(mass**2)', data=d).fit()
m_6_3 = smf.ols('brain ~ mass + I(mass**2) + I(mass**3)', data=d).fit()
m_6_4 = smf.ols('brain ~ mass + I(mass**2) + I(mass**3) + I(mass**4)', data=d).fit()
m_6_5 = smf.ols('brain ~ mass + I(mass**2) + I(mass**3) + I(mass**4) + I(mass**5)', data=d).fit()
m_6_6 = smf.ols('brain ~ mass + I(mass**2) + I(mass**3) + I(mass**4) + I(mass**5) + I(mass**6)', data=d).fit()
m_6_7 = smf.ols('brain ~ 1', data=d).fit()
d_new = d.drop(d.index[-1])
f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(8,3))
ax1.scatter(d.mass, d.brain, alpha=0.8)
ax2.scatter(d.mass, d.brain, alpha=0.8)
for i in range(len(d)):
d_new = d.drop(d.index[-i])
m0 = smf.ols('brain ~ mass', d_new).fit()
# need to calculate regression line
# need to add intercept term explicitly
x = sm.add_constant(d_new.mass) # add constant to new data frame with mass
x_pred = pd.DataFrame({'mass': np.linspace(x.mass.min() - 10, x.mass.max() + 10, 50)}) # create linspace dataframe
x_pred2 = sm.add_constant(x_pred) # add constant to newly created linspace dataframe
y_pred = m0.predict(x_pred2) # calculate predicted values
ax1.plot(x_pred, y_pred, 'gray', alpha=.5)
ax1.set_ylabel('body mass (kg)', fontsize=12);
ax1.set_xlabel('brain volume (cc)', fontsize=12)
ax1.set_title('Underfit model')
# fifth order model
m1 = smf.ols('brain ~ mass + I(mass**2) + I(mass**3) + I(mass**4) + I(mass**5)', data=d_new).fit()
x = sm.add_constant(d_new.mass) # add constant to new data frame with mass
x_pred = pd.DataFrame({'mass': np.linspace(x.mass.min()-10, x.mass.max()+10, 200)}) # create linspace dataframe
x_pred2 = sm.add_constant(x_pred) # add constant to newly created linspace dataframe
y_pred = m1.predict(x_pred2) # calculate predicted values from fitted model
ax2.plot(x_pred, y_pred, 'gray', alpha=.5)
ax2.set_xlim(32,62)
ax2.set_ylim(-250, 2200)
ax2.set_ylabel('body mass (kg)', fontsize=12);
ax2.set_xlabel('brain volume (cc)', fontsize=12)
ax2.set_title('Overfit model')
plt.show()
p = (0.3, 0.7)
-sum(p * np.log(p))
0.6108643020548935
# fit model
m_6_1 = smf.ols('brain ~ mass', data=d).fit()
#compute de deviance by cheating
-2 * m_6_1.llf
94.92498968588757
# standarize the mass before fitting
d['mass_s'] = d['mass'] - np.mean(d['mass'] / np.std(d['mass']))
with pm.Model() as m_6_8 :
a = pm.Normal('a', mu=np.mean(d['brain']), sd=10)
b = pm.Normal('b', mu=0, sd=10)
sigma = pm.Uniform('sigma', 0, np.std(d['brain']) * 10)
mu = pm.Deterministic('mu', a + b * d['mass_s'])
brain = pm.Normal('brain', mu = mu, sd = sigma, observed = d['brain'])
m_6_8 = pm.sample(2000, tune=5000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [sigma, b, a] Sampling 2 chains: 100%|██████████| 14000/14000 [00:10<00:00, 1284.28draws/s]
theta = az.summary(m_6_8)['mean'][:3]
#compute deviance
dev = - 2 * sum(stats.norm.logpdf(d['brain'], loc = theta[0] + theta[1] * d['mass_s'] , scale = theta[2]))
dev
100.60410127191153
This is the original function.
# This function only works with number of parameters >= 2
def sim_train_test(N=20, k=3, rho=[0.15, -0.4], b_sigma=100):
n_dim = 1 + len(rho)
if n_dim < k:
n_dim = k
Rho = np.diag(np.ones(n_dim))
Rho[0, 1:3:1] = rho
i_lower = np.tril_indices(n_dim, -1)
Rho[i_lower] = Rho.T[i_lower]
x_train = stats.multivariate_normal.rvs(cov=Rho, size=N)
x_test = stats.multivariate_normal.rvs(cov=Rho, size=N)
mm_train = np.ones((N,1))
np.concatenate([mm_train, x_train[:, 1:k]], axis=1)
#Using pymc3
with pm.Model() as m_sim:
vec_V = pm.MvNormal('vec_V', mu=0, cov=b_sigma * np.eye(n_dim),
shape=(1, n_dim), testval=np.random.randn(1, n_dim)*.01)
mu = pm.Deterministic('mu', 0 + pm.math.dot(x_train, vec_V.T))
y = pm.Normal('y', mu=mu, sd=1, observed=x_train[:, 0])
with m_sim:
trace_m_sim = pm.sample()
vec = pm.summary(trace_m_sim)['mean'][:n_dim]
vec = np.array([i for i in vec]).reshape(n_dim, -1)
dev_train = - 2 * sum(stats.norm.logpdf(x_train, loc = np.matmul(x_train, vec), scale = 1))
mm_test = np.ones((N,1))
mm_test = np.concatenate([mm_test, x_test[:, 1:k +1]], axis=1)
dev_test = - 2 * sum(stats.norm.logpdf(x_test[:,0], loc = np.matmul(mm_test, vec), scale = 1))
return np.mean(dev_train), np.mean(dev_test)
n = 20
tries = 10
param = 6
r = np.zeros(shape=(param - 1, 4))
train = []
test = []
for j in range(2, param + 1):
print(j)
for i in range(1, tries + 1):
tr, te = sim_train_test(N=n, k=param)
train.append(tr), test.append(te)
r[j -2, :] = np.mean(train), np.std(train, ddof=1), np.mean(test), np.std(test, ddof=1)
2
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 696.53draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 670.08draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 629.84draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 514.13draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 603.98draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 490.54draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 592.44draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 538.58draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 648.23draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 603.95draws/s]
3
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 618.29draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 639.22draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 481.07draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 519.18draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 518.25draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 571.07draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 542.79draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 462.96draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:05<00:00, 390.67draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 684.51draws/s]
4
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 804.55draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 611.91draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 458.24draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 689.49draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 556.63draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 720.40draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 570.87draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 737.06draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 754.59draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 613.29draws/s]
5
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 723.48draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 803.98draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 718.26draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 503.21draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 704.96draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 667.57draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 527.39draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 734.26draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:03<00:00, 548.40draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 667.49draws/s]
6
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 721.66draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 819.60draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:04<00:00, 480.44draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 683.64draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 727.10draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 769.77draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:01<00:00, 1027.99draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 759.98draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 677.01draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [vec_V] Sampling 2 chains: 100%|██████████| 2000/2000 [00:02<00:00, 797.01draws/s]
num_param = np.arange(2, param + 1)
plt.figure(figsize=(10, 6))
plt.scatter(num_param, r[:, 0], color='C0')
plt.xticks(num_param)
for j in range(param - 1):
plt.vlines(num_param[j], r[j,0] - r[j, 1], r[j,0] + r[j,1], color='mediumblue',
zorder=-1, alpha=0.80)
plt.scatter(num_param + 0.1, r[:, 2], facecolors='none', edgecolors='k')
for j in range(param - 1):
plt.vlines(num_param[j] + 0.1, r[j,2] - r[j, 3], r[j,2] + r[j,3], color='k',
zorder=-2, alpha=0.70)
dist = 0.20
plt.text(num_param[1] - dist, r[1, 0] - dist, 'in', color='C0', fontsize=13)
plt.text(num_param[1] + dist, r[1, 2] - dist, 'out', color='k', fontsize=13)
plt.text(num_param[1] + dist, r[1, 2] + r[1,3] - dist, '+1 SD', color='k', fontsize=10)
plt.text(num_param[1] + dist, r[1, 2] - r[1,3] - dist, '+1 SD', color='k', fontsize=10)
plt.xlabel('Number of parameters', fontsize=14)
plt.ylabel('Deviance', fontsize=14)
plt.title(f'N = {n}', fontsize=14)
plt.show()
data = pd.read_csv('Data/cars.csv', sep=',')
with pm.Model() as m_6_15 :
a = pm.Normal('a', mu=0, sd=100)
b = pm.Normal('b', mu=0, sd=10)
sigma = pm.Uniform('sigma', 0, 30)
mu = pm.Deterministic('mu', a + b * data['speed'])
dist = pm.Normal('dist', mu=mu, sd=sigma, observed = data['dist'])
m_6_15 = pm.sample(5000, tune=10000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [sigma, b, a] Sampling 2 chains: 100%|██████████| 30000/30000 [00:31<00:00, 943.55draws/s]
n_samples = 1000
n_cases = data.shape[0]
ll = np.zeros((n_cases, n_samples))
for s in range(0, n_samples):
mu = m_6_15['a'][s] + m_6_15['b'][s] * data['speed']
p_ = stats.norm.logpdf(data['dist'], loc=mu, scale=m_6_15['sigma'][s])
ll[:,s] = p_
n_cases = data.shape[0]
lppd = np.zeros(n_cases)
for a in range(1, n_cases):
lppd[a,] = logsumexp(ll[a,]) - np.log(n_samples)
pWAIC = np.zeros(n_cases)
for i in range(1, n_cases):
pWAIC[i,] = np.var(ll[i,])
- 2 * (sum(lppd) - sum(pWAIC))
412.31568650728855
waic_vec = - 2 * (lppd - pWAIC)
(n_cases * np.var(waic_vec))**0.5
15.162740093813598
d = pd.read_csv('Data/milk.csv', sep=';')
d['neocortex'] = d['neocortex.perc'] / 100
d.dropna(inplace=True)
d.shape
(17, 9)
a_start = d['kcal.per.g'].mean()
sigma_start = d['kcal.per.g'].std()
mass_shared = theano.shared(np.log(d['mass'].values))
neocortex_shared = theano.shared(d['neocortex'].values)
with pm.Model() as m6_11:
alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)
mu = alpha + 0 * neocortex_shared
sigma = pm.HalfCauchy('sigma',beta=10, testval=sigma_start)
kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])
trace_m6_11 = pm.sample(1000, tune=1000)
with pm.Model() as m6_12:
alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)
beta = pm.Normal('beta', mu=0, sd=10)
sigma = pm.HalfCauchy('sigma',beta=10, testval=sigma_start)
mu = alpha + beta * neocortex_shared
kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])
trace_m6_12 = pm.sample(5000, tune=15000)
with pm.Model() as m6_13:
alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)
beta = pm.Normal('beta', mu=0, sd=10)
sigma = pm.HalfCauchy('sigma', beta=10, testval=sigma_start)
mu = alpha + beta * mass_shared
kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])
trace_m6_13 = pm.sample(1000, tune=1000)
with pm.Model() as m6_14:
alpha = pm.Normal('alpha', mu=0, sd=10, testval=a_start)
beta = pm.Normal('beta', mu=0, sd=10, shape=2)
sigma = pm.HalfCauchy('sigma', beta=10, testval=sigma_start)
mu = alpha + beta[0] * mass_shared + beta[1] * neocortex_shared
kcal = pm.Normal('kcal', mu=mu, sd=sigma, observed=d['kcal.per.g'])
trace_m6_14 = pm.sample(5000, tune=15000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [sigma, alpha] Sampling 2 chains: 100%|██████████| 4000/4000 [00:01<00:00, 2441.14draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [sigma, beta, alpha] Sampling 2 chains: 100%|██████████| 40000/40000 [01:46<00:00, 376.87draws/s] There were 14 divergences after tuning. Increase `target_accept` or reparameterize. The acceptance probability does not match the target. It is 0.8907813033282677, but should be close to 0.8. Try to increase the number of tuning steps. The number of effective samples is smaller than 25% for some parameters. Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [sigma, beta, alpha] Sampling 2 chains: 100%|██████████| 4000/4000 [00:02<00:00, 1496.02draws/s] Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [sigma, beta, alpha] Sampling 2 chains: 100%|██████████| 40000/40000 [03:02<00:00, 219.41draws/s] There were 93 divergences after tuning. Increase `target_accept` or reparameterize. There were 23 divergences after tuning. Increase `target_accept` or reparameterize. The number of effective samples is smaller than 25% for some parameters.
az.waic(trace_m6_14, m6_14)
/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:993: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """
waic -16.9226 waic_se 4.87234 p_waic 3.00199 warning True waic_i [-0.5820018909358816, -1.5105480636019253, -1.... waic_scale deviance dtype: object
compare_df = az.compare({'m6_11' : trace_m6_11,
'm6_12' : trace_m6_12,
'm6_13' : trace_m6_13,
'm6_14' : trace_m6_14}, method='pseudo-BMA')
compare_df
/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:993: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:993: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """
waic | p_waic | d_waic | weight | se | dse | warning | waic_scale | |
---|---|---|---|---|---|---|---|---|
m6_14 | -16.9226 | 3.00199 | 0 | 0.957972 | 4.87234 | 0 | True | deviance |
m6_13 | -9.08342 | 1.94869 | 7.83922 | 0.0190146 | 4.05308 | 3.33736 | True | deviance |
m6_11 | -8.72321 | 1.30887 | 8.19943 | 0.0158807 | 3.55127 | 4.60363 | False | deviance |
m6_12 | -7.12231 | 1.99539 | 9.80033 | 0.00713244 | 3.23372 | 4.75175 | False | deviance |
az.plot_compare(compare_df);
diff = np.random.normal(loc=6.7, scale=7.26, size=100000)
sum(diff < 0) / 100000
0.17745
Compare function already checks number of observations to be equal.
coeftab = pd.DataFrame({'m6_11': pm.summary(trace_m6_11)['mean'],
'm6_12': pm.summary(trace_m6_12)['mean'],
'm6_13': pm.summary(trace_m6_13)['mean'],
'm6_14': pm.summary(trace_m6_14)['mean']})
coeftab
m6_11 | m6_12 | m6_13 | m6_14 | |
---|---|---|---|---|
alpha | 0.657407 | 0.364541 | 0.701868 | -1.059950 |
beta | NaN | 0.434416 | -0.030158 | NaN |
beta__0 | NaN | NaN | NaN | -0.095699 |
beta__1 | NaN | NaN | NaN | 2.753594 |
sigma | 0.188607 | 0.191718 | 0.182601 | 0.140068 |
traces = [trace_m6_11, trace_m6_12, trace_m6_13, trace_m6_14]
models = [m6_11, m6_12, m6_13, m6_14]
az.plot_forest(traces, figsize=(10, 5));
kcal_per_g = np.repeat(0, 30) # empty outcome
neocortex = np.linspace(0.5, 0.8, 30) # sequence of neocortex
mass = np.repeat(4.5, 30) # average mass
mass_shared.set_value(np.log(mass))
neocortex_shared.set_value(neocortex)
post_pred = pm.sample_posterior_predictive(trace_m6_14, samples=10000, model=m6_14)
100%|██████████| 10000/10000 [00:19<00:00, 519.32it/s]
milk_ensemble = pm.sample_posterior_predictive_w(traces, 10000,
models,
weights=compare_df.weight.sort_index(ascending=True))
100%|██████████| 10000/10000 [00:19<00:00, 513.64it/s]
plt.figure(figsize=(8, 6))
plt.plot(neocortex, post_pred['kcal'].mean(0), ls='--', color='k')
az.plot_hpd(neocortex, post_pred['kcal'],
fill_kwargs={'alpha': 0},
plot_kwargs={'alpha':1, 'color':'k', 'ls':'--'})
plt.plot(neocortex, milk_ensemble['kcal'].mean(0), color='C1')
az.plot_hpd(neocortex, milk_ensemble['kcal'])
plt.scatter(d['neocortex'], d['kcal.per.g'], facecolor='None', edgecolors='C0')
plt.ylim(0.3, 1)
plt.xlabel('neocortex')
plt.ylabel('kcal.per.g');
/home/osvaldo/anaconda3/lib/python3.7/site-packages/scipy/signal/_arraytools.py:45: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. b = a[a_slice]
import platform
import sys
import IPython
import matplotlib
import scipy
print("""This notebook was created using:\nPython {}\nIPython {}\nPyMC3 {}\nArviZ {}\nNumPy {}\nSciPy {}\nMatplotlib {}\n""".format(sys.version[:5], IPython.__version__, pm.__version__, az.__version__, np.__version__, scipy.__version__, matplotlib.__version__))
This notebook was created using: Python 3.7.1 IPython 6.2.1 PyMC3 3.7.rc1 ArviZ 0.4.0 NumPy 1.15.4 SciPy 1.1.0 Matplotlib 3.0.2