How to generate Multivariate Laplace from Multivariate Normal?

Modeling
1

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);
}
2

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