# Estimate Structured Matrix

**URL:** https://discourse.mc-stan.org/t/estimate-structured-matrix/3353
**Category:** Developers
**Created:** [February 19, 2018, 7:25pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353 "2018-02-19T19:25:11Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![ENFP99](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/enfp99/32/1869_2.png) [@ENFP99](https://discourse.mc-stan.org/u/ENFP99)
#### Post date: [February 19, 2018, 7:25pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353/1 "2018-02-19T19:25:11Z")

</div>

This is the same code, sorry for the duplicate posting, I just keep encountering different types of problems.  
My H matrix, ideally is an indicator matrix filled with 0 and 1, and that is what I really want. When I generate the simulated data, the real H is:  
[,1] [,2] [,3] [,4] [,5] [,6] [,7]  
[1,] 1 1 1 0 1 1 0  
[2,] 1 1 0 0 1 1 0  
[3,] 1 0 1 0 1 1 0  
[4,] 0 1 1 1 1 1 0  
[5,] 1 0 1 0 1 1 0  
[6,] 0 0 0 0 1 1 0  
[7,] 0 0 0 0 1 1 0  
[8,] 1 0 1 0 1 1 0  
[9,] 1 0 0 0 1 1 0  
[10,] 1 1 0 0 1 1 0  
[11,] 1 0 1 1 0 0 1  
[12,] 1 0 1 0 0 0 1  
[13,] 0 1 0 1 0 0 1  
[14,] 1 0 1 0 0 0 1  
[15,] 0 0 0 1 0 0 1  
[16,] 1 1 0 0 0 0 1  
[17,] 0 1 0 0 0 0 1  
[18,] 0 1 0 0 0 0 1  
[19,] 1 0 0 1 0 0 1  
I figured it out now, through asking a series of questions, that Stan cannot estimate such a discrete matrix, then I relaxed my requirement. I used some continuous prior that can introduce sparsity instead of the spike and slab. But as you can see, there is some structure in the H matrix. The first 3 columns have randomly distributed 0 and 1, column 5 and 6 have all 1s for the first few rows and 0s for the bottom few rows and the 7th column has 0s for the first few rows and 1s for the bottom few rows. In order to reflect such a structure on the prior, I used an uniform prior for column 5,6,7. But when I used code like this  
for(j in (k+1):(k+k1))  
for(i in 1:p)  
H[i,j] ~ **uniform(0.9,1)**;  
for(j in (k+1):(k+k1))  
for(i in (p+1):D) H[i,j] ~ **uniform(0,0.1)**;

for(j in (k+k1+1):K)  
for(i in 1:p) H[i,j] ~ uniform(0,0.1);  
for(j in (k+k1+1):K)  
for(i in (p+1):D) H[i,j] ~ uniform(0.9,1);  
It does not work. It will report Log probability evaluates to log(0), i.e. negative infinity."  
" Stan can’t start sampling from this initial value." Only by using  
for(j in (k+1):(k+k1))  
for(i in 1:p)  
H[i,j] ~ uniform(-1,1);  
for(j in (k+1):(k+k1))  
for(i in (p+1):D) H[i,j] ~ **uniform(-1,1);**

for(j in (k+k1+1):K)  
for(i in 1:p) H[i,j] ~ uniform(-1,1);  
for(j in (k+k1+1):K)  
for(i in (p+1):D) H[i,j] ~ **uniform(-1,1)**;  
The code will run. But this is not what I want, the estimate from this set up is far from accurate. I wonder how to fix this problem.  
Attached please is my simulated data and code:

## Generate Data

library(MASS)

p=10; q=9;D=p+q;  
k=4; k1=2; k2=1; K=k+k1+k2;  
n=100;

mu\_z \<- rep(0,K)  
Sigma\_z \<- matrix(0, nrow=K, ncol=K) + diag(K)\*1

Z \<- mvrnorm(n, mu\_z, Sigma\_z)  
H \<- matrix(nrow=D, ncol=K)

set.seed(1)  
for(j in 1:k)  
{H[,j] = rbinom(D, 1, 2/(j+2))}

for(j in (k+1):(k+k1))  
{ for(i in 1:p) H[i,j]=1  
for(i in (p+1):D) H[i,j]=0 }

for(j in (k+k1+1):K)  
{ for(i in 1:p) H[i,j] = 0  
for(i in (p+1):D) H[i,j] = 1}

set.seed(1)  
E\_mux1 \<- runif(K, min =.01, max = 1)  
E\_mux \<- sort(E\_mux1, decreasing = TRUE)  
set.seed(1)  
e\_x \<- runif(K, min = .5, max = 5)  
E\_sigmax \<- diag(K)\*e\_x

set.seed(1)  
E\_x \<- mvrnorm(p, E\_mux, E\_sigmax)  
set.seed(1)  
E\_muy1 \<- runif(K, min =5, max = 40)  
E\_muy \<- sort(E\_muy1, decreasing = TRUE)  
set.seed(1)  
e\_y \<- runif(K, min = .5, max = 20)  
E\_sigmay \<- diag(K)\*e\_y  
set.seed(1)  
E\_y \<- mvrnorm(q, E\_muy, E\_sigmay)  
E \<- rbind(E\_x, E\_y)

Phi \<- E\*H

Theta \<- Z%\*%t(Phi)  
Theta\_x \<- Theta[,1:10]  
Theta\_y \<- Theta[,11:19]

X\_p \<- matrix(nrow=n, ncol=p)  
{for(j in 1:p)  
{ for(i in 1:n)  
{X\_p[i,j] = exp(Theta\_x[i,j])/(1+exp(Theta\_x[i,j]))}}}  
set.seed(1)  
X\<- matrix(nrow=n, ncol=p)  
{for(j in 1:p)  
for(i in 1:n) {X[i,j] = rbinom(1, 2, X\_p[i,j])}}

Y\_mu \<- matrix(nrow=n, ncol=q)  
{for(j in 1:q)  
for(i in 1:n) {Y\_mu[i,j] = Theta\_y[i,j]}}

set.seed(1)  
Y\<- matrix(nrow=n, ncol=q)  
{for(j in 1:q)  
{ for(i in 1:n)  
{Y[i,j] = rnorm(1,Theta\_y[i,j], sd=1/2)}}}

## stan code

library(rstan)  
library(coda)  
SBECCA\_mod \<- '  
data {  
int\<lower=1\> n;  
int\<lower=1\> p;  
int\<lower=1\> q;  
int\<lower=1\> k;  
int\<lower=1\> k1;  
int\<lower=1\> k2;  
int\<lower=1\> D;  
int\<lower=1\> K;  
int\<lower=0, upper=2\> X[n,p];  
real Y[n,q];  
vector[K] mu\_z;  
cov\_matrix[K] Sigma\_z;  
}  
parameters {  
matrix\<lower=-1,upper=1\>[D,K] H;  
matrix[D,K] E;  
matrix[n,K] Z;  
matrix\<lower=0\>[D,k] Pi;  
real\<lower=0\> alpha[K]; # top part of E  
real\<lower=0\> beta[K]; # bottom part of E  
real\<lower=0\> sigma[q]; # modeling data Y, Y is normally distributed  
}

transformed parameters {  
matrix[D,K] Phii;  
matrix[n,D] Theta;  
matrix[n,p] Theta\_x;  
matrix[n,q] Theta\_y;  
matrix[n,p] p\_x;  
matrix[n,q] mu\_y;

Phii \<- E .\* H;  
Theta \<- Z\*transpose(Phii);

for(i in 1:n)  
{for(j in 1:p)  
Theta\_x[i,j] \<- Theta[i,j];

for(l in 1:q)  
Theta\_y[i,l] \<- Theta[i,l+p];}

{for(i in 1:n)  
for(j in 1:p)  
p\_x[i,j] \<- exp(Theta\_x[i,j])/(1+exp(Theta\_x[i,j]));}

{for(i in 1:n)  
for(l in 1:q)  
mu\_y[i,l] \<- Theta\_y[i,l];}  
}  
model {

## Prior

# for Z

for(i in 1:n)  
{Z[i] ~ multi\_normal(mu\_z, Sigma\_z);}

# for alpha, beta (E)

for(j in 1:K)  
{alpha[j] ~ inv\_gamma(0.001, 0.001);  
beta[j] ~ inv\_gamma(0.001, 0.001);} ## This is just an initial value

############################################################################### for H  
for(i in 1:D)  
for(j in 1:k)  
Pi[i,j] ~ exponential(2);  
for(i in 1:D)  
for(j in 1:k)  
H[i,j] ~ normal(0,Pi[i,j]);

for(j in (k+1):(k+k1))  
for(i in 1:p)  
H[i,j] ~ uniform(-1,1);  
for(j in (k+1):(k+k1))  
for(i in (p+1):D) H[i,j] ~ uniform(-1,1);

for(j in (k+k1+1):K)  
for(i in 1:p) H[i,j] ~ uniform(-1,1);  
for(j in (k+k1+1):K)  
for(i in (p+1):D) H[i,j] ~ uniform(-1,1);  
############################################################################ for H

for(j in 1:K)  
for(i in 1:p)  
E[i,j] ~ normal(0, alpha[j]);

for(j in 1:K)  
for(l in 1:q)  
E[l,j] ~ normal(0, beta[j]);

# for variance of y

for(j in 1:q)  
sigma[j] ~ inv\_gamma(0.01, 0.01);

## Generate Data

for(i in 1:n)  
{  
for(j in 1:p)  
{  
X[i,j] ~ binomial(2,p\_x[i,j]);  
}  
for(l in 1:q)  
{  
Y[i,l] ~ normal(mu\_y[i,l], sigma[l]);  
}  
}  
}  
'  
SBECCA\_dat\<-list(“X”=X,“Y”=Y, “mu\_z” = rep(0,7), “Sigma\_z” = diag(rep(1,7)),“n”=100,“p”=10,“q”=9,“k”=4, “k1”=2, “k2”=1, “D”= 19, “K”=7)  
happyfit3 \<- stan(model\_code = SBECCA\_mod, data = SBECCA\_dat, thin=3, iter = 100, chains =2)

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [February 20, 2018, 7:22pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353/2 "2018-02-20T19:22:07Z")

</div>

If you have uniforms, the parameters need to be declared with the same bounds. Otherwise, unconstrained values will lead to illegal constrained values, and Stan will freeze because it can’t find parameter values randomly that satisfy constraints or it’ll revert to random walks when it steps out of support.

You can use three back-ticks to

```
```
```
 set off code. One before and one after.

Stan’s not really set up to do this kind of hard choice. Aki’s been writing parameters on softer forms of variable selection.

---

<div class="post-metadata">

### Author: ![ENFP99](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/enfp99/32/1869_2.png) [@ENFP99](https://discourse.mc-stan.org/u/ENFP99)
#### Post date: [February 20, 2018, 7:30pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353/4 "2018-02-20T19:30:33Z")

</div>

Can you elaborate this part a little bit??

You can use three back-ticks to

```
set off code. One before and one after.
```

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [February 20, 2018, 9:37pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353/5 "2018-02-20T21:37:45Z")

</div>

You got the quote before the code, but not the one after. You want it to look like this:

```
```
double a = 10;
a *= 3;
```
```

and then it’ll render as

```
double a = 10;
a *= 3;

```

---

<div class="post-metadata">

### Author: ![ENFP99](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/enfp99/32/1869_2.png) [@ENFP99](https://discourse.mc-stan.org/u/ENFP99)
#### Post date: [February 22, 2018, 6:28pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353/6 "2018-02-22T18:28:46Z")

</div>

Do you think it is possible to marginalize H matrix element-wise given that H has such a structure?

I received a lot of suggestions on marginalizing discrete parameters on the previous question I asked, but it is a simplified situation, the real situation I need to deal with is like this (structured H). Do you think it is possible to marginalize this?

I know that this is technically not a Stan question, but I just feel you guys have a much broader knowledge than merely “stan”, would like to get some suggestion from you! .^ ^.

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [February 23, 2018, 8:05pm UTC](https://discourse.mc-stan.org/t/estimate-structured-matrix/3353/7 "2018-02-23T20:05:44Z")

</div>

When the discrete parameters are organized along with the data, then it’s possible. When the discrete parameters interact with each (i.e., they’re not all conditionally independent), as in variable selection, they can’t be marginalized out efficiently enough to be coded in Stan.