//Bayesian nonparametric finite population inference when the survey outcome is binary
//reference: Si, Y., Pillai, N. and Gelman, A. (2015), Bayesian nonparametric weighted sampling inference, Bayesian Analysis, 10(3), 605-625.
//Author: Yajuan Si
//Latest edit date: 03-25-2017
data {                //input data elements 
  real<lower=1> N;       //population size, define it as real to be used in n/N
  int<lower=1> n;       //sample size
  int<lower=1> J;       //number of poststratification cells
  vector<lower=1>[J] vec_n;      //vector of sample cell sizes
  vector[J] vec_w;      //vector of cell weight values
  int vec_y[J];          //cell total of survey bonary response
}
transformed data {
  vector[J] x;          //covariates in GP regression, logarithm of the weights
  vector[J] vec2_y;     //change vec_y to a vector, convenience for vector operation
  x <- log(vec_w); 
  for (j in 1:J){
    vec2_y[j] <- vec_y[j];
  }
}
}
parameters {
  real beta_raw;               //reparameterization of beta: beta=beta_raw * sd_beta (2.5)
  real<lower=0> eta;    //partial sill in GP covariance kernel
  real<lower=0> kappa;       //decay parameter in GP covariance kernel
  real<lower=0, upper=1> delta_w; //weight assigned to independent prior specification
  vector<lower=0>[J] vec_N_p; //population cell size proportiong N_j/N
  vector[J] z; //linear combination predictor in the logistric regression with the coronical parameter
}
transformed parameters {
  real beta;  //coefficients in the mean function of GP
  real eta_sq; //square of eta
  real c;      //normalizing constant 
  for (j in 1:J)
  beta <- 2.5 * beta_raw;
  eta_sq <-  pow(eta, 2); 
  c <-  sum(vec_N_p) / N * n / sum(vec_N_p ./ vec_w); 
}
model {
// 1) y|w under GP prior
  vector[J] mu0; //mean of mu
  matrix[J,J] Sigma; //covariance matrix
  mu0 <- x * beta; 
 for (i in 1:J) {
    for (j in i:J) {
      Sigma[i,j] <- eta_sq * (1-delta_w) * exp(-kappa * fabs(x[i] - x[j])); //squared exponential kernel
      Sigma[j,i] <- Sigma[i,j];
    }
    Sigma[i,i] <- Sigma[i,i] + eta_sq * delta_w; //combined with independent prior 
  }
  z ~ multi_normal(mu0, Sigma);
  vec_y ~ binomial_logit(vec_n, z); //logistic regression
  // 2) n_j|w  
  vec_n ~ poisson(c * (vec_N_p ./ vec_w)); 

  //this can be changes as multinomial distribution when J is small
  n ~ poisson(c * sum(vec_N_p ./ vec_w));
  vec_n ~ multinomial(vec_N_p ./ vec_w/sum(vec_N_p ./ vec_w));

  // 3) prior for hyperparameters, weakly informative
  eta ~ cauchy(0,5);
  kappa ~ gamma(0.5,0.5); 
  beta_raw ~ normal(0,1); 
}
generated quantities {//predictions
  vector[J] theta; //probability
  real y_est;
  y_est <- 0;
  for (j in 1:J){
    theta[j] <- inv_logit(z[j]);
  }
  y_est <- y_prd - sum(theta .* vec_n- vec2_y) / N /sum(vec_N_p) ;
}