# How to generate Multivariate Laplace from Multivariate Normal?

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 Zero = rep_vector(0, 3);
}
parameters {
cholesky_factor_corr Lcorr;
vector<lower=0> sigma;
row_vector alpha[J];  // player skills
}
transformed parameters {
// define correlation and covariance matrices
corr_matrix R;
cov_matrix 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);
}

The responsibility of the model block is to take a set of parameters and calculate the log-prior and log-likelihood. You can think of this as just a function: given parameters, calculate a value. Mathematically, for parameter \theta and data y, Bayes rule can be expressed as

p(\theta \mid y) \propto \underset{\text{prior}}{\underbrace{p(\theta)}} \times \underset{\text{likelihood}}{\underbrace{p(y \mid \theta)}}

All the model block does is compute the log of the right hand side: given \theta, compute \log p(\theta) + \log p(y \mid \theta). Note that this is a deterministic function, there’s no randomness involved here – that’s why you can’t use any of the random number generating functions in the transformed parameters or model blocks.

Another important note is that the following syntax

alpha ~ multi_normal(Zero, Sigma);


Is just shorthand for

target += multi_normal_lpdf(alpha | Zero, Sigma)


The full syntax might make it more clear that there are no random numbers being generated here – we’re just incrementing the calculated log-prior plus log-likelihood by the log-probability density function of the multivariate normal distribution evaluated at alpha with parameters Zero and Sigma.

To get to the specifics of your question, if I understand correctly you would like to assign a multivariate Laplace prior to alpha. To do this, you will need to be able to compute the log-probability density function of the multivariate Laplace distribution. I don’t think this is built into Stan (like the log-density of the multivariate Normal distribution is), so you will need to code it yourself.

The basic idea is that you write a new function (say multi_laplace_lpdf) that calculates the log-density of the multivariate Laplace distribution. Then, in your model block, you would be able to just write

alpha ~ multi_laplace(Zero, Sigma);


More details and an example can be found in the Stan manual: 20.5 User-defined probability functions | Stan User’s Guide.

Hope that helps!

1 Like