Confirm that drawing parameters in Python matches those in R

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)









    



Using tmb_tmp/linreg.cpp.
Compiled in 24.2s.

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()









    



iter: 1  value: 9.618947 mgc: 327.9179 ustep: 1 
iter: 2  mgc: 2.797762e-14 
iter: 1  mgc: 2.797762e-14 
Matching hessian patterns... Done
outer mgc:  7.140877 
iter: 1  mgc: 2.104983e-13 
iter: 1  mgc: 2.104983e-13 
outer mgc:  1.651893 
iter: 1  mgc: 3.703704e-13 
iter: 1  mgc: 2.566836e-13 
iter: 1  mgc: 2.566836e-13 
outer mgc:  0.1627763 
iter: 1  mgc: 2.566836e-13 
iter: 1  mgc: 2.566836e-13 
outer mgc:  0.01984795 
iter: 1  mgc: 2.580158e-13 
iter: 1  mgc: 2.580158e-13 
outer mgc:  0.0002043639 
iter: 1  mgc: 2.553513e-13 
iter: 1  mgc: 2.553513e-13 
outer mgc:  2.530986e-07 
iter: 1  mgc: 2.4869e-13 
iter: 1  mgc: 2.4869e-13 
outer mgc:  3.235856e-12 
iter: 1  mgc: 2.4869e-13 

Model optimization complete in 0.1s.


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

iter: 1  mgc: 2.4869e-13 
outer mgc:  3.235856e-12 
iter: 1  mgc: 2.633449e-13 
outer mgc:  0.01598401 
iter: 1  mgc: 2.606804e-13 
outer mgc:  0.01601601 

Simulated 100 draws in 0.2s.

alpha:
	mean	[ 0.93484804]
	sd	[ 0.19260968]
	draws	[[ 0.84137591  1.06381489 ...,  0.97978633  0.96196277]]
	shape	(1, 100)
Beta:
	mean	[ 1.00298141]
	sd	[ 0.0360791]
	draws	[[ 1.03785315  0.96909651 ...,  0.98399896  1.01262166]]
	shape	(1, 100)

Simulate parameters

Imports



In [9]:

    
import scipy.sparse.linalg
import scikits.sparse.cholmod
from scipy.sparse import csc_matrix
import numpy as np
from PyMB import get_R_attr
r = m.R.r

Means $\mu$



In [10]:

    
R_mu = get_R_attr(m.sdreport, 'par.random')
py_mu = np.array(R_mu)

Joint precision $\Sigma^{-1}$



In [11]:

    
joint_prec_full = get_R_attr(m.sdreport, 'jointPrecision')
R_joint_prec = r('function(mat, ran) { ii <- rownames(mat) %in% ran; return(mat[ii,ii]) }')(joint_prec_full, m.random)
py_joint_prec = np.array(r['as.matrix'](R_joint_prec))

$z$



In [12]:

    
R_z = r['matrix'](r['rnorm'](100*py_mu.shape[0]), ncol=100)
py_z = np.array(R_z)

$\mathrm{Chol}(\Sigma^{-1})$



In [14]:

    
R_chol = r['Cholesky'](R_joint_prec)
py_chol = scikits.sparse.cholmod.cholesky(csc_matrix(py_joint_prec))

$L^{\prime}z$

Note: would typically use scikits.sparse.cholmod.cholesky.solve_Lt, but there seems to be a bug there: https://github.com/njsmith/scikits-sparse/issues/9#issuecomment-76862652



In [16]:

    
R_Lz = r['solve'](r['t'](r['as'](R_chol, 'Matrix')), R_z)
py_Lz = scipy.sparse.linalg.spsolve(py_chol.L().T, py_z)

$PL^{\prime}z$



In [17]:

    
R_PLz = r['solve'](r['as'](R_chol, 'pMatrix'), R_Lz)
py_PLz = py_chol.apply_Pt(py_Lz)

Draws



In [18]:

    
R_draws = r('function(mu,PLz) { return(mu + PLz) }')(R_mu, R_PLz)
py_draws = py_mu.reshape([py_mu.shape[0],1]) + py_PLz

Confirm equality



In [21]:

    
np.allclose(np.array(r['as.matrix'](R_draws)), py_draws)









    Out[21]:





True