Define Dirichlet parameters in advance settings

Modeling
1

Hello,

I am trying to fit a model with Dirichlet prior. However, my parameter specification is quite complex. I have three parameters, let’s say w, r and m. These three follows Dirichlet prior. For all these three parameters should be defined i times. In-other words, my parameters are wi, ri and mi where i=1:400. I can fit this prior easily as follows. (theta_1 = (wi, ri, mi))

Data{
row_vector [3] beta1;
}
parameters {
    
     simplex [3] theta_1 [400];
  
}

model {
  // Define priors  
  theta_1 ~ dirichlet(beta1);
  
dat ~ tenismodel(theta_1);

}

I have done this. Now I want an advance setting than above. I want to fit parameters, theta_1 = (wi ^j, ri^j, mi^j) . Now I have both superscript j and subscript i. As an example, when i = 1 and j = 2 I want to have parameters:

[w1^1, w1^2, r1^1, r1^2, m1^1, m1^2]

I tried to explain my problem as much as I can. If anyone have an idea how we can implement these parameter settings in stan please let me know. I really appreciate your help.

Thank You.

2

I’m not 100% sure from your description what you need. Is it this (for an illustrative case where j ranges from 1 to 10): simplex[3] theta_1[400, 10]?

Note that in newer versions of Stan, this would be written as array[400, 10] simplex[3] theta_1

3

Thanks jsocolar for the reply. Does this simplex[3] theta_1[400, 10] means, each parameter repeats 10 times for all 400 cases?.

4

This means that there are 4000 3-simplexes, with each simplex indexed by i and j

5

Thanks, I will run this and check is that what I want.

6

Hi Jsocolar,

I tried this method and still I am getting an error. “Ill-typed arguments to ‘~’ statement. No distribution ‘dirichlet’ was found with the correct signature.”

Think this is because my beta1. I supplied it as an vector in R. b1 <- c(rep(0.0078, 81), rep(0.0913,81), rep(0.9009,81)). How do we need to supply beta 1 to the model?. Is that a matrix?.

7

It’s hard to say for sure since you are just showing parts of your model definition, but most straightforwardly, beta1 would also be a 400 x 81 array of 3-vectors.

8

Hi Jsocolar,

This following is my full code.


// Create user define function

functions{

  real tenismodel_lpdf(matrix dat, vector[] theta, vector[] alpha, int[] ServerID, int[] ReceiverID, int[] t1,
  int[] t2, int[] t3){

  vector[rows(dat)] prob;
  real  t;
  real  x;
  int s;
  int re;
  real out;
  real a;
  real f;
  real ws;
  real rs;
  real ms;
  real wr;
  real rr;
  real mr;
  real pr;
    
  for (i in 1:rows(dat)){
    
    t <- dat[i, 1];
    x <- dat[i, 2];
    s <- ServerID[i];
    re <- ReceiverID[i];
    a <- theta[s, 1];
    f <- theta[s, 2];

    
    if(x == -1){
      
      prob[i] <- f; 
      
    }else if(t == 1){
      
      prob[i] <- a^x * f^(1 - x);

    }else if(fmod(t, 2) == 0){
     
     vector [t1[i]] p1 ;

     for(j in 1:t1[i]) {
       
        rs <- alpha[s, j];
        rr <- alpha[re, j];
        
        p1[j] <- (rr)^(2*j) * (rs)^(2*j + 1);
     }
     
       pr <-  prod(p1);
      
        ws <- alpha[s, 1];
        ms <- alpha[s, 3];
        wr <- alpha[re, 1];
        mr <- alpha[re, 3];
      prob[i] <- ((1- a - f) * pr * (mr)^x * (wr)^(1-x)); 
    
    }else {
      
      vector [t2[i]] p1;
     
     for(j in 1:t2[i]) {
       
        rs <- alpha[s, j];
        rr <- alpha[re, j];
        
        p1[j] <- (rr)^(2*j);
     }
     
     vector [t3[i]] p2;
     
     for(j in 1:t3[i]) {
       
        rs <- alpha[s, j];
        rr <- alpha[re, j];
        
        p2[j] <- (rs)^(2*j + 1);
     }
     
     pr <- prod(p1) * prod(p2);
     ws <- alpha[s, 1];
     ms <- alpha[s, 3];
     wr <- alpha[re, 1];
     mr <- alpha[re, 3];
     
       prob[i] <- ((1- a - f) * pr * (ms)^(1 - x) * (ws)^x);
    }
  }

  out <- sum(log(prob));
  return out;

 }
}


// The input data is a vector 'y' of length 'N'.
data {
  int N; // number of observations (rally lengths)
  int M; // number of players
  int t1; // number of touches
  int t2; // number of touches
  int t3; // number of touches
  matrix [N, 2] dat; // touches, indicator, servereID and receiverID 
  int ServerID [N];
  int ReceiverID [N];
  int <lower = 2> m; 
  row_vector [3] beta1;
  row_vector [3] beta2;


}

parameters {
    
     simplex [m] theta_1 [M, 81];
     
     simplex [m] theta_2 [M, 81];
  
}

// The model to be estimated.
model {
  // Define priors

  
  theta_1 ~ dirichlet(beta1);

  theta_2 ~ dirichlet(beta2);

  dat ~ tenismodel(theta_1, theta_2, ServerID, ReceiverID, t1, t2, t3);

}

To feed values for the model, I use the following R code:


fit_model3 <- stan(file = 'R/model_3.stan',
                   data = list(N = len, M = num_players,
                               dat = as.matrix(men_data),
                               ServerID = ServerID,
                               ReceiverID = ReceiverID,
                               m = 3,
                               beta1 = b1,
                               beta2 = b2,
                               t1 = t1,
                               t2 = t2, t3 = t3), 
                   chains = 1, iter = 1000)

where my b1 and b2 are:

vector_1 <- rep(0.0078, 81)
vector_2 <- rep(0.0913,81)
vector_3 <- rep(0.9009,81)
b1<- array(c(vector_1,vector_2, vector_3))
b2<- array(c(vector_1,vector_2, vector_3))

I do not think the way I defined b1 and b2 is correct.

9

Which dirichlet parameters are meant to apply to which elements of theta_1?

10

For theta_1, m is number of parameters, which is 3; M are the number of players and 81 is how many time each parameter repeats. According to my example, i = M and j = 81.

11

Right now theta_1 involves m\times M \times 81 parameters. Are you now saying that it should contain only 3 parameters, or do you still mean 3 parameters for each of the m \times M elements? Do you wish for all of those M \times m simplexes to all have the same Dirichlet prior whose parameters are given by the single 3-vector beta_1?

12

Yes, I have three parameters in my model. As an example I have three parameters a, b, c. I want a to have aij where i = M and j = 81. Also, bij again i = M and j = 81 and cij where i = M and j = 81. Is this possible? . I think this statement is true ‘’ Do you wish for all of those
M
×
m
simplexes to all have the same Dirichlet prior whose parameters are given by the single 3-vector beta_1?‘’

13

I think the most straight forward solution here is not to worry about vectorizing and to write

for(i in 1:m){
  for(j in 1:M){
    theta_1[i,j] ~ dirichlet(beta_1);
  }
}
14

Hi Jscolar,


data {
  int N; // number of observations (rally lengths)
  int M; // number of players
  int t1; // number of touches
  int t2; // number of touches
  int t3; // number of touches
  int K;
  matrix [N, 2] dat; // touches, indicator, servereID and receiverID 
  int ServerID [N];
  int ReceiverID [N];
  int <lower = 2> m; 
  row_vector [3] beta1;
  row_vector [3] beta2;


}

parameters {
    
     simplex [m] theta_1 [M];
     
     simplex [m] theta_2 [M];
  
}

// The model to be estimated.
model {
  // Define priors

  for(i in 1:M){
  for(j in 1:K){
    theta_1[i,j] ~ dirichlet(beta1);
  }
}

  for(i in 1:M){
  for(j in 1:K){
    theta_2[i,j] ~ dirichlet(beta2);
  }
}


  dat ~ tenismodel(theta_1, theta_2, ServerID, ReceiverID, t1, t2, t3);

}

No I gave beta1 and beta2 as vectors in R. “beta1 ← c(0.0815, 0.3814, 0.5371)
beta2 ← c(0.0715, 0.7424, 0.1861)”.

But I am getting the same error. The error is:

Ill-typed arguments to '~' statement. No distribution 'dirichlet' was found with the correct signature.

What am I doing in correct when defining Dirichlet parameters?. My M is number of players; which is i and my K is j ; how many times each parameter repeats. That means each of the three parameters in beta1 is repeating K times.

15

You have a few problems here.

On the error, the signature for the dirichlet distribution takes vectors, not row vectors.

Once you resolve the above compile error, and assuming the rest is ok to compile, you will also have a runtime error for dimension mismatch between your data declaration row_vector[3] beta1 and what you are giving the model with b1, which is a longer than 3 vector.

16

From the doc:

The Dirichlet probability functions are overloaded to allow the simplex θ 
and prior counts (plus one) α to be vectors or row vectors (or to mix the 
two types).
1 Like
17

Ah, I think you have theta as an array of simplexes,

but you loop through theta twice, leaving only a real for the dirichlet, which expects the simplex,

for(i in 1:M){
  for(j in 1:K){
    theta_1[i,j] ~ dirichlet(beta1);
  }
}
1 Like
18

@niro previously you had theta as a 2-dimensional array of 3-simplexes, but now you have a 1-dimensional array of 3-simplexes. Which do you intend?

19

Jscolar, I think I should explain my question again. May be it is not clear here. This is what I want.

Let’s say I have four players in my data set, i.e M = 4, and I need to estimate three parameters for each of these players (m = 3). Let’s say they are a, b, c. These three parameters should repeat for each player 2 times (k = 2). So my “a” would be;

player 1 - a11, a12
player2- a21, a22
player3- a31, a32
player4- a41, a42

Also, for next parameter b;

player 1 - b11, b12
player2- b21, b22
player3- b31, b32
player4- b41, b42

For the last one c also as above.

So I need to set prior distribution of these to Dirichlet distribution. This is my problem.

20

Yup, so in your most recent Stan code, you declared a one-dimensional array of simplexes, but what you want is a two-dimensional array of simplexes.