Raising a multi_normal function (bivariate normal) model to a power

I’m very new to Stan so I’m sorry if this is a basic question, but I want to run a model that’s simply a bivariate normal distribution raised to power but I can’t figure out how to make it work. So something basically like:


Here’s the full code I’ve tried to run (I’ve tried a few subtle variations to this as well

// The input data is a vector 'y' of length 'N'.
data {
  int<lower=0> N;
  vector[2] y[N];

// The parameters accepted by the model. 
parameters {
  vector[2] mu; //means
  vector<lower=0>[2] sigma; //sds
  real<lower = -1,upper = 1> rho; //correltation
  real<lower = 0,upper = 1> alpha;

transformed parameters{
  cov_matrix[2] Sigma; // create first correlation matrix
  Sigma[1,1] = square(sigma[1]); // variance 
  Sigma[2,2] = square(sigma[2]); // varance
  Sigma[1,2] = rho*sigma[1]*sigma[2]; // covariance
  Sigma[2,1] = Sigma[1,2]; // covariance

// The model to be estimated. We model the output
// 'y' to bebivariate normally distributed with mean 'mu'
// and  covariance  'Sigma'.
model {
mu[1] ~ normal(-4,1);
mu[2] ~ normal(25,3);
sigma[1] ~ cauchy(1,.2);
sigma[2] ~ cauchy(3,.2);
rho ~ normal(-.2,.1);
alpha ~ normal(.5,.1);

This code works perfectly until I try to raise multi_normal(mu, Sigma) to a power.


I can’t tell for sure what you want the power function to achieve.

Do you want to encode the assumption that y arises from a probability distribution with PDF proportional to normalPDF^alpha?

Or do you want to encode that y arises from a probability distribution proportional to exp(normal_lpdf^alpha)?

Or do you mean that y^(-alpha) arises from a normal distribution? (note that this is a completely different thing from the other two options; you cannot treat the sampling statement ~ like an equals sign and maintain the same model as you perform identical operations on both sides of the ~.

So I just want alpha to scale the existing bivariate probability density function. For example, y ~ multi_normal(mu,Sigma)^.5 should be more spread out than y~multi_normal(mu,Sigma)^1. Does that make sense? In r I can just do something like sum(dmvnorm(x,mu,Sigma))^.5. I think I don’t fully understand how the multi_normal functions work in Stan though.

Note that the product of Gaussian PDFs is proportional to another Gaussian, so you might be able to do this all analytically and then just sample from a new normal distribution. However, barring that, the most straightforward way to achieve this is with target += alpha * multi_normal_lpdf(y | mu, Sigma).

1 Like

that worked, thanks so much!