# Poor sample result and slow sampling speed on Multivariate normal distribution

**URL:** <https://discourse.pymc.io/t/poor-sample-result-and-slow-sampling-speed-on-multivariate-normal-distribution/9207>\
**Category:** v3\
**Tags:** modeling\
**Created:** [April 14, 2022, 1:48pm UTC](https://discourse.pymc.io/t/poor-sample-result-and-slow-sampling-speed-on-multivariate-normal-distribution/9207 "2022-04-14T13:48:10Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![qipengchen](https://avatars.discourse-cdn.com/v4/letter/q/65b543/32.png) [@qipengchen](https://discourse.pymc.io/u/qipengchen)\
**Post date:** [April 14, 2022, 1:48pm UTC](https://discourse.pymc.io/t/poor-sample-result-and-slow-sampling-speed-on-multivariate-normal-distribution/9207/1 "2022-04-14T13:48:10Z")

</div>

Hi! I’m estimating a longitudinal item response theory model, and the speed of sampling is very slow.  
My model is as follows:

\theta\_{tn} \sim MN(\mu\_{\theta}, \Sigma\_\theta)

logit(P\_{tni})= \alpha\_i\theta\_{tn}+\beta\_{i} 

Y\_{tni} \sim Bernoulli(P\_{tni})

Y\_{tni} is the response of student n(1,…,N) on item i(1,…,I),time point t(1,…,T).

For identify of the model, I need to constrain the mean and the variance of first time point. So I’m constrain the first element of the mean vector to be 0, and the first element of the covariance matrix to be 1. like that:

\mu\_\theta = (0,\mu\_{\theta2},...,\mu\_{\theta T})

And disassemble the covariance matrix as diagonal standard deviation matrix S and correlation matrix \Omega ,then constrain the first element of diagonal standard deviation matrix as 1.

\mathbf\Sigma = \textbf S \mathbf\Omega \textbf S

\textbf S = \left(\begin{array}{l} 1\\ & \sigma\_{\theta 2}\\ & & \ddots& \\ & && \sigma\_{\theta T} \end{array}\right)

 \mathbf\Omega= \left(\begin{array}{l} 1 &\rho\_{12}& \cdots & \rho\_{1T}\\ \rho\_{12} &1 &\cdots &\rho\_{2T} \\ \vdots &\vdots &\ddots &\vdots\\ \rho\_{1T}& \rho\_{2T}&\cdots & 1 \end{array}\right)\\ 

The observation variable Y is generated by the following code

```auto
T = 5
N = 300
I = 10
theta_mean = np.zeros(T)
theta_cov = [[1. , 0.9 , 0.81 , 0.729 , 0.6561],
            [0.9 , 1. , 0.9 , 0.81 , 0.729],
            [0.81 , 0.9 , 1. , 0.9 , 0.81],
            [0.729 , 0.81 , 0.9 , 1. , 0.9],
            [0.6561, 0.729 , 0.81 , 0.9 , 1.]]
theta = np.random.multivariate_normal(theta_mean,theta_cov,size=(N))
beta = np.random.normal(0,1,(1,1,I))
alpha = np.random.normal(1,0.15,(1,1,I))
logit_p = alpha*theta.T.reshape(T,N,1)+beta
Y = np.random.binomial(1,np.exp(logit_p)/(1+np.exp(logit_p)))

```

Pymc3 model code. LKJcorr as the prior of correlation matrix.

```auto
triu_idx = np.triu_indices(T,k=1)
with pm.Model() as long_irt:
    corr_array = pm.LKJCorr("corr", n=T, eta=1.0)
    cor_tri = tt.set_subtensor(tt.zeros((T,T))[triu_idx],corr_array)
    corr_mat = cor_tri + tt.eye(5) + cor_tri.T 
    # The standard deviation of theta, constrain sigma[0,0]=1
    sigma = tt.concatenate([[1],pm.HalfNormal("s",1,shape=(T-1))]) 
    sigma_mat = tt.diag(sigma)
    cov = pm.Deterministic("cov",(sigma_mat).dot(corr_mat).dot(sigma_mat))
    # mean of MvNormal, constrain theta_mu[0]=0
    theta_mu = tt.concatenate([[0],pm.Normal("theta_mu",0,10,shape=(T-1))])
    
    theta = pm.MvNormal("theta",theta_mu,cov=cov,shape=(N,T))
    
    alpha_mu = pm.HalfNormal("alpha_mu",5)
    beta_mu = pm.Normal("beta_mu",0,5)
    beta_sigma = pm.HalfNormal("beta_sigma",5)
    alpha_sigma = pm.HalfNormal("alpha_sigma",5)
    alpha = pm.Lognormal("alpha",alpha_mu,alpha_sigma,shape=(1,1,I))
    beta = pm.Normal("beta",beta_mu,beta_sigma,shape=(1,1,I))
    
    pm.Bernoulli("Y",logit_p=theta.T.reshape((T,N,1))*alpha+beta,observed=Y)

```

With this pymc3 model, I got poor sample result and slow sampling speed.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/d/d1ba7978a0577494e32946ccf1b1389066ae84ae.png)

But, the model seems work well when I setting the independent normal distribution as the priors of theta. I also noticed that the tutorials setting prior like this.[Fitting a Multidimensional Item-Response Theory Model in Stan](https://rfarouni.github.io/assets/projects/BayesianMIRT/BayesianMIRT.html)

```auto
with pm.Model() as long_irt:
    theta = pm.Normal("theta",0,1,shape=(N,T))

    alpha_mu = pm.HalfNormal("alpha_mu",5)
    beta_mu = pm.Normal("beta_mu",0,5)
    beta_sigma = pm.HalfNormal("beta_sigma",5)
    alpha_sigma = pm.HalfNormal("alpha_sigma",5)
    alpha = pm.Lognormal("alpha",alpha_mu,alpha_sigma,shape=(1,1,I))
    beta = pm.Normal("beta",beta_mu,beta_sigma,shape=(1,1,I))
    
    pm.Bernoulli("Y",logit_p=theta.T.reshape((T,N,1))*alpha+beta,observed=Y)

```

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/f/f5a5c8b3a6482383a0018dd51b986bffb670ae1a.png)

However, I need to modeling correlation of the theta by covariance matrix structure. **Why MvNormal sampling so slow and how can I get valid samples from multivariate normal distribution in this model?**

---

<div class="post-metadata">

**Author:** ![chartl](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/chartl/32/1515_2.png) [@chartl](https://discourse.pymc.io/u/chartl)\
**Post date:** [April 14, 2022, 6:50pm UTC](https://discourse.pymc.io/t/poor-sample-result-and-slow-sampling-speed-on-multivariate-normal-distribution/9207/2 "2022-04-14T18:50:17Z")

</div>

I have a strong suspicion that this all stems from the use of a centered parameterization. What happens if instead of using a MvNormal you did something like:

```auto
L_mat = cor_tri + tt.diag(sigma)
theta_offsets = pm.normal('theta_offsets', 0, 1, shape=(T, N))
theta = pm.Deterministic('theta', theta_mu + tt.dot(L_mat, theta_offsets))

```

which will have row-covariance LL^T = \Sigma and row means \theta\_\mu.

---

<div class="post-metadata">

**Author:** ![qipengchen](https://avatars.discourse-cdn.com/v4/letter/q/65b543/32.png) [@qipengchen](https://discourse.pymc.io/u/qipengchen)\
**Post date:** [April 15, 2022, 12:39pm UTC](https://discourse.pymc.io/t/poor-sample-result-and-slow-sampling-speed-on-multivariate-normal-distribution/9207/3 "2022-04-15T12:39:39Z")

</div>

Thanks, It helps me a lot! I never thought that multivariate distribution could use non-centered parameterization.  
I would like to know more about non-centered parameterization, is there any recommended reference material?

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [April 16, 2022, 3:26am UTC](https://discourse.pymc.io/t/poor-sample-result-and-slow-sampling-speed-on-multivariate-normal-distribution/9207/4 "2022-04-16T03:26:53Z")

</div>

> [@qipengchen](#):
>
> I would like to know more about non-centered parameterization, is there any recommended reference material?

The STAN user guide includes [an excellent discussion on reparameterization](https://mc-stan.org/docs/2_29/stan-users-guide/reparameterization.html), including centered/non-centered. [This blog post](https://betanalpha.github.io/assets/case_studies/identifiability) by Michael Betancourt might also interest you, it touches on reparameterization in a broader discussion of model identification, sampler degeneracy, and what to do when NUTS is going slow or diverging.

For PyMC specific, the [case study on multilevel modeling example notebook](https://docs.pymc.io/en/v3/pymc-examples/examples/case_studies/multilevel_modeling.html) has examples of centered and non-centered implementations for the univariate and multivariate cases. I got a lot out of this notebook, but it’s a bit hidden!