What's wrong with my stan code for estimating parameters with multivariate distribution

Modeling
1

Hello, everyone. I have a problem that X,Y \sim N(\mu,\Sigma), where the TRUE value is that \mu=(10,20) and \Sigma=(\begin{matrix}10&4\\4&30\end{matrix}).

I KNOW how to estimate with stan-defined functions and now i want to practice in user-defined functions. Now here is my stan code:

code='
  functions{
    real fun_log(vector X,vector Y,real mu1,real mu2,real sigma1, real sigma2, real rho){
    vector [num_elements(X)] prob;
    real lprob;
    for (i in 1:num_elements(X)){
     ##Likelihood function
      prob[i] <- (2*pi()*sigma1*sigma2*sqrt(1-rho^2))^(-1)*exp(-1/2/(1-rho^2)*((X[i]-mu1)^2/sigma1^-2*rho*(X[i]-mu1)*(Y[i]-mu2)/sigma1/sigma2+(Y[i]-mu2)^2/sigma2^2));
    }
    lprob <- sum(log(prob));
    return lprob;
    }
  }
data{
  int N;
  vector[N] X;
  vector[N] Y;
}
parameters{
  real<lower=0> rho;
  real<lower=0> sigma1;
  real<lower=0> sigma2;
  real mu1;
  real mu2;
}
model{
append_row(X,Y)~fun(rho,sigma1,sigma2,mu1,mu2);
}
'

But it says that

Unknown variable: append_row
No matches for: 
  vector ~ fun(real, real, real, real, real)

Can some one help me to debug the model? Thanks!

2

Judging from the error message, maybe you need to declare a local variable for the return value of append_row, and have this as the target of your sampling statement?

Also, your return from the function is a single real, so it won’t accept a vector as target?

3

I’m pretty sure that the problem is that fun_log is defined incorrectly. In a log-probability function, all arguments except the first are the parameters to the function. Here, Y is treated as a parameter, so it belongs on the right-hand side of the “~” operator, as part of the arguments to fun. In short, given how you defined fun_log, the corresponding sampling statement has to be

X ~ fun(Y, rho,sigma1,sigma2,mu1,mu2);

What you probably really want to do is define fun_log like this:

real fun_log(matrix XY, real mu1,real mu2,real sigma1, real sigma2, real rho){
    vector [cols(XY)] prob;
    real lprob;
    for (i in 1:cols(XY)){
     ##Likelihood function
      prob[i] <- (2*pi()*sigma1*sigma2*sqrt(1-rho^2))^(-1)*exp(-1/2/(1-rho^2)*((XY[1,i]-mu1)^2/sigma1^-2*rho*(XY[1,i]-mu1)*(XY[2,i]-mu2)/sigma1/sigma2+(XY[2,i]-mu2)^2/sigma2^2));
    }
    lprob <- sum(log(prob));
    return lprob;
    }

Then define the data block as follows:

data{
  int N;
  matrix[2, N] XY;
}

The sampling statement in your model would then be

XY ~ fun(rho,sigma1,sigma2,mu1,mu2);
1 Like