My tennis model below involves modeling player skills using a multivariate normal distribution, where there are 3 components to a player’s skill: his skill on hard, clay, and grass courts. An extension I’d like to try is modeling their skill as a multivariate Laplace distribution.

One way to transform a multivariate normal sample to a multivariate Laplace sample is the following: if Y \sim N_3(0, \Sigma) and W \sim Exp(1), then \sqrt{W}Y=X is multivariate Laplace of dimension 3 with mean 0 and covariance matrix \Sigma. This is noted at the bottom of the multivariate Laplace wikipedia page here.

My problem is that I cannot use Stan’s `exponential_rng`

function in either the `model`

or `transformed parameters`

block. I tried doing it `generated quantities`

which you can see below, but I don’t think it’s correct, as the `generated quantities`

block is only evaluated **after** a posterior sample has been drawn, meaning the player skills `alpha`

are not multivariate normal anymore.

Does anyone have any ideas?

```
data {
int<lower=0> J; // total number of players
int<lower=0> N_hard; // number of matches on hard surface
int<lower=0> N_clay; // on clay surface
int<lower=0> N_grass; // on grass surface
matrix[N_hard, J] X_skill_hard; // design matrix for the skills on hard
matrix[N_clay, 2*J] X_skill_clay; // design matrix for the skills on clay
matrix[N_grass, 2*J] X_skill_grass; // design matrix for the skills on grass
int winner_y_hard[N_hard]; // number of sets won by the winner for match n on hard
int winner_y_clay[N_clay]; // number of sets won by the winner for match n on clay
int winner_y_grass[N_grass]; // number of sets won by the winner for match n on grass
int total_n_hard[N_hard]; // number of total sets played in match n on hard
int total_n_clay[N_clay]; // number of total sets played in match n on clay
int total_n_grass[N_grass]; // number of total sets played in match n on grass
}
transformed data {
vector[3] Zero = rep_vector(0, 3);
}
parameters {
cholesky_factor_corr[3] Lcorr;
vector<lower=0>[3] sigma;
row_vector[3] alpha[J]; // player skills
}
transformed parameters {
// define correlation and covariance matrices
corr_matrix[3] R;
cov_matrix[3] Sigma;
R = multiply_lower_tri_self_transpose(Lcorr); // does Lcorr %*% t(Lcorr)
Sigma = quad_form_diag(R, sigma); // does sigmaI %*% R %*% sigmaI
vector[3*J] alpha_vector = to_vector(to_matrix(alpha));
// Is there a way to generate a multivariate Laplace sample from my alpha's here?
}
model {
// priors
Lcorr ~ lkj_corr_cholesky(1);
sigma ~ cauchy(0, 4);
alpha ~ multi_normal(Zero, Sigma);
// Or could I transform the alpha's here?
// sample
winner_y_hard ~ binomial_logit(total_n_hard, X_skill_hard*(alpha_vector[1:J]));
winner_y_clay ~ binomial_logit(total_n_clay, X_skill_clay*(append_row(alpha_vector[1:J], alpha_vector[(J+1):(2*J)])));
winner_y_grass ~ binomial_logit(total_n_grass, X_skill_grass*(append_row(alpha_vector[1:J], alpha_vector[(2*J+1):(3*J)])));
}
generated quantities {
row_vector[J] rate = rep_row_vector(1, J); // (1, 1, 1, ..... , 1) J-vector of rates, one for each players skill vector
vector[J] W = to_vector(exponential_rng(rate));
matrix[J, 3] alpha_laplace = rep_matrix(sqrt(W), 3) .* to_matrix(alpha);
vector[3*J] alpha_laplace_vector = to_vector(alpha_laplace);
}
```