In [1]:

    

%pylab inline
import numpy as np
from scipy import stats
import pymc3 as pm
import seaborn as sns
import matplotlib.pyplot as plt


    

    




Populating the interactive namespace from numpy and matplotlib






    




WARNING (theano.configdefaults): Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

In [2]:

    

data = np.random.randn(100)


    

In [3]:

    

with pm.Model() as model: 
    mu = pm.Normal('mu', mu=0, sd=1, testval=0)
    sd = pm.HalfNormal('sd', sd=1)
    n = pm.Normal('n', mu=mu, sd=sd, observed=data)


    

    




Applied log-transform to sd and added transformed sd_log_ to model.

In [4]:

    

means, sds, elbos = pm.variational.advi(model=model, n=20000, accurate_elbo=True)


    

    




Iteration 0 [0%]: ELBO = -586.43
Iteration 2000 [10%]: Average ELBO = -386.82
Iteration 4000 [20%]: Average ELBO = -157.07
Iteration 6000 [30%]: Average ELBO = -136.69
Iteration 8000 [40%]: Average ELBO = -131.44
Iteration 10000 [50%]: Average ELBO = -130.67
Iteration 12000 [60%]: Average ELBO = -130.66
Iteration 14000 [70%]: Average ELBO = -130.67
Iteration 16000 [80%]: Average ELBO = -130.66
Iteration 18000 [90%]: Average ELBO = -130.66
Finished [100%]: Average ELBO = -130.66

In [5]:

    

with model:
    step = pm.NUTS()
    trace = pm.sample(1000, step)


    

    




100%|██████████| 1000/1000 [00:00<00:00, 2698.98it/s]

In [6]:

    

print(trace['mu'].mean())
print(trace['sd_log_'].mean())
print(trace['mu'].std())
print(trace['sd_log_'].std())


    

    




0.00766113915397
-0.146682044436
0.0897388737967
0.0773784687122

In [7]:

    

ax = sns.distplot(trace['mu'], label='NUTS')
xlim = ax.get_xlim()
x = np.linspace(xlim[0], xlim[1], 100)
y = stats.norm(means['mu'], sds['mu']).pdf(x)
ax.plot(x, y, label='ADVI')
ax.set_title('mu')
ax.legend(loc=0)


    

    




/usr/local/var/pyenv/versions/miniconda3-4.0.5/envs/testenv/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[7]:





<matplotlib.legend.Legend at 0x14339a7b8>






    

In [8]:

    

ax = sns.distplot(trace['sd_log_'], label='NUTS')
xlim = ax.get_xlim()
x = np.linspace(xlim[0], xlim[1], 100)
y = stats.norm(means['sd_log_'], sds['sd_log_']).pdf(x)
ax.plot(x, y, label='ADVI')
ax.set_title('log(sigma)')
ax.legend(loc=0)


    

    




/usr/local/var/pyenv/versions/miniconda3-4.0.5/envs/testenv/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[8]:





<matplotlib.legend.Legend at 0x1434bb278>






    

In [9]:

    

import matplotlib.pyplot as plt
plt.plot(elbos[5000:])


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0x1437d35f8>]






    

In [ ]: