//Bayesian nonparametric finite population inference when the survey outcome is continuous //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 int N; //population size int n; //sample size int J; //number of poststratification cells vector[J] vec_n; //vector of sample cell sizes vector[J] vec_w; //vector of cell weight values vector[J] vec_y; //vector of cell means for survey response: \bar{y}_j vector[J] vec_s2; //vector of cell sum of squared error for survey response: s_j^2 } transformed data { vector[J] x; //covariates in GP regression, logarithm of the weights x <- log(vec_w); } parameters { real beta_raw; //reparameterization of beta: beta=beta_raw * sd_beta (2.5) real sigma_sq; //nugget in GP covariance kernel real eta; //partial sill in GP covariance kernel real kappa; //decay parameter in GP covariance kernel real delta_w; //weight assigned to independent prior specification vector[J] vec_N_p; //population cell size proportiong N_j/N vector[J] mu; //realizations of mean functions } 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 vector[J] sigma_j; 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 } mu ~ multi_normal(mu0, Sigma); //GP realizations //---use cell mean and sum of squared errors as sufficient statistics for normal distributions---// sum(vec_s2) / sigma_sq ~ chi_square(n-1); lp__ <- lp__ - log(sigma_sq); //Jacobian transformation sigma_j <- sqrt(sigma_sq ./ vec_n); vec_y ~ normal(mu, sigma_j); //---above is used for better computation performance---// // 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 beta_raw ~ normal(0,1); eta ~ cauchy(0,5); kappa ~ gamma(0.5,0.5); sigma_sq ~ inv_gamma(1,0.5); } generated quantities { //predictions real y_est; y_est <- y_prd - sum((mu - vec_y) .* vec_n) / N / sum(vec_N_p); }