Stochastic Blockmodel by Stan

Modeling
1

I am trying to sample from the density given by the last expression on page 5.
Here is my Stan code. It doesn’t converge.

(upload://8IPHO5hhKDQ5E1Hk00jUPVudszd.pdf) (544.2 KB) ](https://arxiv.org/pdf/1203.3083.pdf)

data{
  int N;
  int K;
  int X[N*N];
  
}

parameters{
  vector<lower=0>[K] alpha1; //各オブジェクトの混合割合
  vector<lower=0>[K] alpha2;
  simplex[K] pi1; //K面サイコロ(Dirichret分布)
  simplex[K] pi2;
  real<lower=0> a0; 
  real<lower=0> b0;
  vector<lower=0,upper=1>[K*K] theta; //クラスタ間の関係の強さ
}

model{
  real lp[K*K];

  pi1 ~ dirichlet(alpha1);
  pi2 ~ dirichlet(alpha2);
  theta ~ beta(a0, b0);

  for (k in 1:K){
    for (l in 1:K){
      lp[k+ K*(l-1)] = log(pi1[k]) + log(pi2[l]) + bernoulli_lpmf(X | theta[k+ K*(l-1)]);
    }
  }
  target += lp;
}

I’m new to Stan. I would like you to help coding this by Stan properly.

1 Like
2

The weirdest part to me in the model formulation is bernoulli_lpmf(X | theta[k+ K*(l-1)]) are you sure you want to express the whole X as a mixture of components, each individually bernoulli? This would mean that the data can be fully summarized by sum(X) and N (your code would be equivalent to binomial_lpmf(sum(X) | N * N, theta[k+ K*(l-1)]))?

In its current form, the model likely suffers from non-identifiability ie. that there are multiple model configurations that lead to the same (or very similar) mode (maximum) of the log density, because if you swap pi1 with pi2 and transpose theta , you’ll get exactly the same likelihood.

A few additional notes:

  • It is useful to provide more details about your issues (the exact messages, the data you pass to the model)
  • Stan has types for matrices and 2D arrays, so you could have matrix<lower=0,upper=1>[K,K] theta or int X[N,N] making it easire to write your model
  • Useful information can be extracted diagnostic plots from ShinyStan or as described in the Bayesplot Vignette https://cran.r-project.org/web/packages/bayesplot/vignettes/visual-mcmc-diagnostics.html

Best of luck!