Bivariate Alpha-Skew Normal model

Modeling
1

Greetings folks

I am a new user of the Stan programing language, and I could use some help with the implementation of a bivariate response variable following the Alpha-Skew Normal distribution.

For instance, let’s consider (Y_1,Y_2) \sim \text{bi-ASN}([\mu_1,\mu_2]',\Sigma,[\alpha_1,\alpha_2]') and the pdf is

f(y_1,y_2|\mu_1,\mu_2,\sigma_1,\sigma_2,\rho,\alpha_1,\alpha_2)=\frac{1 + (1-\alpha_1(\frac{y_1-\mu_1}{\sigma_1})-\alpha_2(\frac{y_2-\mu_2}{\sigma_2}))^2}{K}\phi_{\mathbb{Y}}(y_1,y_2)

whereas K=2+\alpha_1^2+\alpha_2^2+2\alpha_1\alpha_2\rho and the \phi_{\mathbb{Y}}(\cdot) is the bivariate normal dist. density (further details, Louzada, Ara & Fernandes, 2016).

functions{
  real biASN_log(matrix y, vector mu, vector a, real s1, real s2, real r){
    real prob = 0;
    for (i in 1:(num_elements(y))){
      prob += log((1-a[1]*((y[i][1]-mu[1])/s1)-a[2]*((y[i][2]-mu[2])/s2))^2+1)-
      log(s1*s2*sqrt(1-r^2)*(2+a[1]^2+a[2]^2+2*a[1]*a[2]*r))-
      ((((y[i,1]-mu[1])/s1)^2+((y[i,2]-mu[2])/s2)^2-2*r*((y[i,1]-mu[1])/s1)*((y[i,2]-mu[2])/s2))/(1-r^2));
    }
    return prob;
  }
}

data { 
  int<lower=0> N;    // number of data items
  matrix[N,2] Y;     // outcome matrix
}

parameters {
  row_vector[2] mu;  // Location parameter
  row_vector[2] a;   // Asymmetry parameter (alpha)
  real<lower=0> s1;  // Sqrt of variances for Y1
  real<lower=0> s2;  // Sqrt of variances for Y2
  real<lower=0> r;   // correlation parameter
} 

model {
  // Priors
  to_vector(mu) ~ normal(0,10);
  to_vector(a) ~ exponential(0.5);
  // Likelihood
  Y ~ biASN(mu,a,s1,s2,r);
}

However, I am getting an error in the model’s likelihood line “Available argument signatures for biASN: Real return type required for probability function”.

Thanks in advance for the support!

1 Like
2

Welcome to the Discourse. I believe that you can only use the sampling statements (~) with functions ending with _lpdf and _lpmf. So try changing the function name to biASN_lpdf. Alternatively, you could keep the name as is and increment target.

target += biASN_log(Y, mu,a,s1,s2,r);

See 7.4 Sampling statements | Stan Reference Manual

1 Like
3

Thanks for the reply. The tip mentioned above was a great point, which is now fixed. Nevertheless, I believe that dimensionally operations are not enabling the compilation. The code is now like:

// Bivariate ASN model
functions{
  real biASN_log(matrix Y, vector mu0, matrix sigma0, vector alpha0){
    vector[rows(Y)] prob;
    real lprob;
    for (i in 1:rows(Y)){
      prob[i] = log( (1-alpha0[1]*((Y[i,1]-mu0[1])/sigma0[1,1])-alpha0[2]*((Y[i,2]-mu0[2])/sigma0[2,2]))^2 +1 ) 
      - log(2+alpha0[1]^2+alpha0[2]^2+2*alpha0[1]*alpha0[2]*sigma0[1,2]) + multi_normal_lpdf(Y[i,:] | mu0, sigma0);
    }
    lprob = sum(prob);
    return lprob;
  }
}

data { 
  int<lower=0> n;
  vector[2] Y[n];
}
parameters {
  vector[2] mu;
  vector[2] alpha;
  vector<lower=0>[2] lambda;
  real<lower=-1,upper=1> r;
} 
transformed parameters {
  vector<lower=0>[2] sigma;
  cov_matrix[2] T;

  // Reparameterization
  sigma[1] = inv_sqrt(lambda[1]);
  sigma[2] = inv_sqrt(lambda[2]);
  T[1,1] = square(sigma[1]);
  T[1,2] = r * sigma[1] * sigma[2];
  T[2,1] = r * sigma[1] * sigma[2];
  T[2,2] = square(sigma[2]);
}
model {
  // Data
  Y ~ biASN(mu, T, alpha);
}
4

Note that you changed Y from a matrix to an array of vectors. But your function still requires a matrix as the first argument. So if you change Y back to a matrix (and replace _log with _lpdf), you should be good to go.

A few side-notes: Stan no longer accepts the old array specification that you’re using (vector[2] Y[n]). The new version looks like this: array[n] vector[2] Y.

Second, you might get a tiny speedup by specifying T as a 2-by-2 matrix rather than a covariance matrix (e.g. matrix[2,2 T). Your covariance matrix should be valid by design, so you don’t need the additional check that Stan does to confirm. That probably isn’t too relevant here, but it may be helpful later on as you build more complex models. Removing checks that you know are valid can be helpful.

Finally, I assume that you are specifying the covariance matrix with lambda and r for a reason (e.g. adding priors to those components later). But, as-is, you could just specify the covariance matrix in the parameters block and not have to worry about recreating it yourself. Of course, in that case, you’d want to use the cov_matrix data structure to ensure that it is positive definite.

1 Like
5

Thank you very much. It worked! =)

1 Like