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.