PyMB Example

Import Module



In [1]:

    
import sys
sys.path.append('../..')
import PyMB

Create a new model



In [2]:

    
m = PyMB.model(name='linreg')

Define the model

$$ Y \sim \mathcal{N}(\hat{Y},\sigma) $$$$ \hat{Y} = \alpha + B x $$

Based on https://github.com/kaskr/adcomp/blob/master/tmb_examples/linreg.cpp

See the TMB tutorial for more information on writing custom models



In [3]:

    
linreg_code = '''
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() (){
// DATA
  DATA_VECTOR(Y);
  DATA_VECTOR(x);

// PARAMETERS
  PARAMETER(alpha);
  PARAMETER(Beta);
  PARAMETER(logSigma);

// MODEL
  vector<Type> Y_hat = alpha + Beta*x;
  REPORT(Y_hat);
  Type nll = -sum(dnorm(Y, Y_hat, exp(logSigma), true));
  return nll;
}
'''

Compile

See additional compiler options such as setting the R and TMB library directories via help(PyMB.model.compile)



In [4]:

    
m.compile(codestr=linreg_code)









    



Saving model to tmb_tmp/linreg.cpp.
Compiled in 24.3s.

Simulate data



In [5]:

    
import numpy as np
m.data['x'] = np.arange(10)
m.data['Y'] = m.data['x'] + 0.5 + np.random.rand(10)

Set initial parameter values



In [6]:

    
m.init['alpha'] = 0.
m.init['Beta'] = 0.
m.init['logSigma'] = 0.

Set random parameters

The model likelihood will be integrated wrt the random parameters. See here for more information.



In [7]:

    
m.random = ['alpha', 'Beta']

Fit the model



In [8]:

    
m.optimize()









    



Matching hessian patterns... Done

Model optimization complete in 0.1s.


--------------------------------------------------------------------------------


Simulated 100 draws in 0.2s.

alpha:
	mean	[ 1.20647766]
	sd	[ 0.12282757]
	draws	[[ 1.12340849  1.42198335 ...,  1.25513521  1.38185193]]
	shape	(1, 100)
Beta:
	mean	[ 0.94143064]
	sd	[ 0.02300771]
	draws	[[ 0.96776205  0.90909791 ...,  0.92343494  0.89284158]]
	shape	(1, 100)

Extract fitted values



In [9]:

    
print(m.report('Y_hat'))









    



[ 1.20647766  2.14790831  3.08933895  4.0307696   4.97220024  5.91363089
  6.85506153  7.79649218  8.73792282  9.67935346]

Examine joint density



In [10]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
df = pd.DataFrame({ k: v['draws'][0] for k,v in m.parameters.iteritems() })
g = sns.PairGrid(df, diag_sharey=False)
g.map_lower(sns.kdeplot, cmap='Blues_d')
g.map_upper(plt.scatter)
g.map_diag(sns.kdeplot, lw=3)









    



/usr/bin/anaconda/lib/python2.7/site-packages/matplotlib/axes/_axes.py:476: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "