Defining dirichlet distributions for segments of a matrix

Modeling
, ,
1

I have a cell-mean model where I am trying to put a Dirichlet prior on the design matrix.
I am really struggling to implement this in RStan and would really appreciate some help.
I am struggling to just get a basic inflexible implementation to work.

Here is my current model without the Dirichlet with the desired model (hard coded/inflexible) for the example data provided below (commented out due to invalid syntax):

data {
  int<lower=0> N;   // number of data items
  int<lower=0> K;   // number of means
  matrix[N,K] x;   // design matrix
  vector[N] y;      // outcome vector
}
parameters {
  vector[K] mu;       // coefficients for predictors
  real<lower=0> sigma;  // error scale
  real eta[K];        // Error in mean
  //simplex[5] x_w1; // Weights in mean estimation
  //simplex[5] x_w2; // Weights in mean estimation
}
//transformed parameters{
//  matrix[N,K] x_w;
//  x_w[1:5,1] = x_w1; x_w[6:10,1] = 0;
//  x_w[1:5,2] = 0; x_w[6:10,1] = x_w2;
//} Should have been set according to x
model {
//  x_w1 ~ dirichlet(1/2); Uniformative priors for simplicty
//  x_w2 ~ dirichlet(1/2);
// 
  y ~ normal(x * mu, sigma); //Should have been x_w instead of x
  sigma ~ gamma(1, 1);
  for (i in 1:K){
     mu  ~ normal(sigma*eta[i], sigma);
  }
  eta ~ normal(0, 1);
}

Simulated data:

# Simulate data in R
set.seed(1)
y <- rnorm(10)
y[6:10] <- y[6:10] + rnorm(5,1, seq(.5, 1.5, length.out = 5))
x <- model.matrix(~0+factor(rep(c(1,2), each = 5)))
dat <- list(
    N = 10,
    K = 2,
    y = y,
    x = x
)

Any help or suggestions (alternative model designs) are highly appreciated.

  • Windows
  • R package version 2.21.5
2

Sorry for the slow response. I don’t fully understand what the intended relationship between x and x_w is (or is x supposed to be completely removed from the working model?). Can you show the code using x_w that isn’t working and explain what happens when you try to use it?

3

Thank you for replying. The intention is to produce a weighted mean estimator (using inferred priors on the dirichlet) for each mean. x just gives the positions that corresponds to values for a particular mean.
E.g if x is the following:

x <- model.matrix(~0+factor(rep(c(1,2), each = 3)))
colnames(x) <- paste0('c', 1:2)
x
  c1 c2
1  1  0
2  1  0
3  1  0
4  0  1
5  0  1
6  0  1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$`factor(rep(c(1, 2), each = 3))`
[1] "contr.treatment"

Then the indexes for x_{1i}=1 should be replaced by the dirichelt and similar for x_{2i}.
Let me know if this brought clarity or if you still like some pseudo ish code!

4

I’m still not sure I understand; tell me if this is right.

You want a 10 x 2 matrix X_w such that X_w[1:5, 1] \sim Dirichlet(...), X_w[6:10, 2] \sim Dirichlet(...), and all other elements are zero?

To achieve this, declare two parameter vectors with the desired Dirichlet priors and then stick them into a matrix of zeros in transformed parameters. But I’m not certain this is what you intend.

1 Like