Matrix of Parameters: best way to specify

Hi all,

I’m using Stan for the first time and I’d like to get feedback on the way I specified my model. I’m trying to predict whether two proteins interact given a particular vector-space encoding. Each pair of features needs a parameter informing whether that pair of features increases or decreases the probability of interaction. I was able to specify the model and run it, but I’m unhappy with the need for the for-loops to specify the prior and likelihood. Is there a way around this? The stan code and PyStan code are below. Thanks!

data {
  int<lower=0> P; // number of proteins
  int<lower=0> N; // number of features
  matrix[P, N] feature_vectors; // data matrix
  int<lower=0,upper=1> interaction_matrix[P, P]; // 0-1 matrix indicating interaction/lack of interaction
parameters {
  matrix[N,N] w; // matrix of feature-feature weights contributing to interaction
transformed parameters {
  real interaction_strength[P, P];
  for (p1 in 1:P)
    for (p2 in 1:P)
      interaction_strength[p1, p2] = feature_vectors[p1] * w * feature_vectors[p2]';
model {
  for (p1 in 1:P)
    for (p2 in 1:P)
      interaction_matrix[p1, p2] ~ bernoulli_logit(interaction_strength[p1, p2]);
  for (n1 in 1:N)
    for (n2 in 1:N)
      w[n1,n2] ~ normal(0, 1);
import pystan
import numpy as np

model_fn = 'interaction_model.stan'

feature_vectors = np.array([
    [.3, 0, 0, .9],
    [.2, .5, .5, 0],
    [.8, .2, 0, 0]
interaction_matrix = np.array([
    [0, 0, 0],
    [0, 0, 1],
    [0, 1, 0]
P, N = feature_vectors.shape
dat = {
    'P': P,
    'N': N,
    'feature_vectors': feature_vectors,
    'interaction_matrix': interaction_matrix
sm = pystan.StanModel(model_fn)
fit = sm.sampling(data=dat, iter=1000, chains=4)

For w, you can to

to_vector(w) ~ normal(0,1);

The other one, you can just change the kernel inside the outer loop

for (p1 in 1:P) {
  vector[P] eta;
  for (p2 in 1:P) eta[p2] = interaction_strength[p1, p2];
  interaction_matrix[p1, ] ~ bernoulli_logit(eta);

Thanks! The new code for w looks much nicer. Is your suggestion for interaction_matrix compiled in a way that makes it faster or is there some other reason to use it? It looks equally verbose to the original.

Faster. Verboseness in a compiled language is not that big a deal as long as the code is clear.