skew.probito_cod <- '
  data {
    int<lower=0> k; 			// // predict
    int<lower=0> n; 			// observat numbers (number of documents)
    int<lower=0,upper=1> y[n]; 	       // response variable
    matrix[n,k] X;  			
    }
transformed data { 
  vector[rows(csr_extract_w(X))] wX;
  int vX[size(csr_extract_v(X))];
  int uX[size(csr_extract_u(X))]; 
  // generate sparse matrix representation of X
  wX = csr_extract_w(X); 
  vX = csr_extract_v(X); 
  uX = csr_extract_u(X);
} 
parameters {
    real beta0;
    vector[k] beta_std;			
    real<lower=-1,upper=1> delta;
    real<lower=0> esc1_global;
    real<lower=0> esc2_global;
    vector<lower=0>[k] esc_local;
  }
 transformed parameters {
  vector[k] beta;
  real alfa;
  real<lower=0> tau; 
  vector<lower=0>[k] lambda;
  vector[n] prob;
  vector[n] eta; 
  vector[n] mu; 
  alfa = delta/sqrt(1-pow(delta,2));
  tau = esc1_global* sqrt(esc2_global);
  lambda = esc_local;
  beta= beta_std .* sqrt(lambda)*tau;
  eta= csr_matrix_times_vector(rows(X), cols(X), wX, vX, uX, beta);  
  mu= beta0 + eta;
  for(i in 1:n){
  prob[i] = skew_normal_cdf(mu[i],0,1,alfa);
 } 
  }
  model {    
    esc1_global ~ normal(0, 1);
    esc2_global ~ inv_gamma(0.5, 0.5);
    esc_local ~ exponential(1);
    beta0 ~ normal(0, 10);
    beta_std ~ normal(0, 1);
    delta ~ uniform(-1, 1);
    y ~ bernoulli(prob);
}
generated quantities { 
vector[n] log_lik; 
for (i in 1:n){
log_lik[i] = bernoulli_lpmf(y[i] | prob[i]);
}
}
'