How to model a multivariate multilevel linear model

Questions
1

Here is my code, model in a native way, Is there a way to rewrite this code with less variable? Thanks…

with pm.Model() as model_4:

  grp_alpha_mu = pm.Normal('grp_alpha_mu', mu=0, sd=10)
  grp_alpha_sd = pm.HalfCauchy('grp_alpha_sd', 5)

  grp_beta_mu_0 = pm.Normal('grp_beta_mu_0', mu=0, sd=10)
  grp_beta_sd_0 = pm.HalfCauchy('grp_beta_sd_0', 5)
  grp_beta_mu_1 = pm.Normal('grp_beta_mu_1', mu=0, sd=10)
  grp_beta_sd_1 = pm.HalfCauchy('grp_beta_sd_1', 5)
  grp_beta_mu_2 = pm.Normal('grp_beta_mu_2', mu=0, sd=10)
  grp_beta_sd_2 = pm.HalfCauchy('grp_beta_sd_2', 5)
  grp_beta_mu_3 = pm.Normal('grp_beta_mu_3', mu=0, sd=10)
  grp_beta_sd_3 = pm.HalfCauchy('grp_beta_sd_3', 5)
  grp_beta_mu_4 = pm.Normal('grp_beta_mu_4', mu=0, sd=10)
  grp_beta_sd_4 = pm.HalfCauchy('grp_beta_sd_4', 5)
  grp_beta_mu_5 = pm.Normal('grp_beta_mu_5', mu=0, sd=10)
  grp_beta_sd_5 = pm.HalfCauchy('grp_beta_sd_5', 5)

  alpha = pm.Normal('alpha', mu=grp_alpha_mu, sd=grp_alpha_sd, shape=n_site)

  beta_0 = pm.Normal('beta', mu=grp_beta_mu_0, sd=grp_beta_sd_0, shape=n_site)
  beta_1 = pm.Normal('beta_1', mu=grp_beta_mu_1, sd=grp_beta_sd_1, shape=n_site)
  beta_2 = pm.Normal('beta_2', mu=grp_beta_mu_2, sd=grp_beta_sd_2, shape=n_site)
  beta_3 = pm.Normal('beta_3', mu=grp_beta_mu_3, sd=grp_beta_sd_3, shape=n_site)
  beta_4 = pm.Normal('beta_4', mu=grp_beta_mu_4, sd=grp_beta_sd_4, shape=n_site)
  beta_5 = pm.Normal('beta_5', mu=grp_beta_mu_5, sd=grp_beta_sd_5, shape=n_site)
  sigma = pm.HalfCauchy('sigma', 5)

   mu = pm.Deterministic('mu', alpha[site_idx] + beta_0[site_idx]*X_train_siteResort['辐照度'] + beta_1[site_idx]*X_train_siteResort['风速'] +
                     beta_2[site_idx]*X_train_siteResort['风向'] + beta_3[site_idx]*X_train_siteResort['温度'] +
                     beta_4[site_idx]*X_train_siteResort['湿度'] + beta_5[site_idx]*X_train_siteResort['压强'])
  y_pred = pm.Normal('y_pred', mu=mu, sd=sigma, observed=y_train_siteResort.values)
1 Like
2

This is a little bit cleaner. You can subsequently use a for loop to build up mu:

import theano
import theano.tensor as tt

with pm.Model() as tight2:
    beta_mu_hp = pm.Normal('beta_mu_hp', 0., 1., shape=7)
    beta_sd_hp = pm.HalfCauchy('beta_sd_hp', 1., shape=7)
    
    alpha = pm.Normal('alpha', mu=beta_mu_hp[0], sd=beta_sd_hp[0], shape=n_site)
    
    beta_mat = [pm.Normal('_beta_{}'.format(i), mu=beta_mu_hp[i+1], sd=beta_sd_hp[i+1], shape=n_site) 
                for i in range(6)]
    beta = pm.Deterministic('beta', tt.stack(beta_mat, axis=1))
    
    sigma = pm.HalfCauchy('sigma', 1)
    
    tr_tight = pm.sample_prior_predictive(100)

@junpenglao

I feel like this should work:

with pm.Model() as tight:
    beta_mu_hp = pm.Normal('beta_mu_hp', 0., 1., shape=7)
    beta_sd_hp = pm.HalfCauchy('beta_sd_hp', 1., shape=7)
    
    alpha = pm.Normal('alpha', mu=beta_mu_hp[0], sd=beta_sd_hp[0], shape=n_site)
    
    beta = pm.Normal('beta', mu=beta_mu_hp[1:].reshape((6,1)), sd=beta_sd_hp[1:].reshape((6,1)), shape=(6,n_site))
    sigma = pm.HalfCauchy('sigma', 1)
    
    tr_tight = pm.sample_prior_predictive(100)

but it generates an error during prior predictive sampling:

~/anaconda3/lib/python3.7/site-packages/pymc3-3.6-py3.7.egg/pymc3/distributions/distribution.py in <listcomp>(.0)
    670     elif broadcast_shape[:len(size_tup)] == size_tup:
    671         suffix = broadcast_shape[len(size_tup):] + dist_shape
--> 672         samples = [generator(*args, **kwargs).reshape(size_tup + (1,)) for _ in range(np.prod(suffix, dtype=int))]
    673         samples = np.hstack(samples).reshape(size_tup + suffix)
    674     else:

ValueError: cannot reshape array of size 600 into shape (100,1)
3

@lucianopaz it seems to be a shape issue for prior predicted - would your PR fix this?

4

Ouch! It looks like a dumb bug in generate_samples. I’ll try to fix it as soon as I can.

2 Likes
5

Thanks chartl, I’am pretty new to the pymc and theano, how to using this model to my data to get the
trace? now I write the mu like this:

mu = pm.Deterministic('mu', alpha[site_idx] + beta_mat[0][site_idx]*all_site_train['辐照度'] + beta_mat[1][site_idx]*all_site_train['风速'] +
                     beta_mat[2][site_idx]*all_site_train['风向'] + beta_mat[3][site_idx]*all_site_train['温度'] +
                     beta_mat[4][site_idx]*all_site_train['湿度'] + beta_mat[5][site_idx]*all_site_train['压强'])
y_pred = pm.Normal('y_pred', mu=mu, sd=sigma, observed=y_train_scaled)

still a naive way…

6

For the trace it’s just pm.sample().

Also you can clan up the code a little bit if you have something like

sitenames = ['辐照度', '风速', ...]
beta_eff = 0
for beta_idx, site_name in enumerate(sitenames):
  beta_eff += beta_mat[beta_idx][site_idx]*all_site_train[site_name]
mu = pm.Deterministic('mu', alpha[site_idx] + beta_eff)
7

thanks chartl, It’s really helpful