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

Modeling
,
1

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:

y~multi_normal(mu,Sigma)^alpha;

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 {
y~multi_normal(mu,Sigma)^alpha;
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.

Thanks!

2

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 ~.

3

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.

4

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
5

that worked, thanks so much!