# Defining dirichlet distributions for segments of a matrix

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 x_w1; // Weights in mean estimation
//simplex 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

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?

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
attr(,"contrasts")
attr(,"contrasts")\$`factor(rep(c(1, 2), each = 3))`
 "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!

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