Implement Gibbs 2PNO Sampler
Gibbs_2PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, burnin, chain_length = 10000L)
| Y | A N by J |
|---|---|
| mu_xi | A two dimensional |
| Sigma_xi_inv | A two dimensional identity |
| mu_theta | The prior mean for theta. |
| Sigma_theta_inv | The prior inverse variance for theta. |
| burnin | The number of MCMC samples to discard. |
| chain_length | The number of MCMC samples. |
Samples from posterior.
# simulate small 2PNO dataset to demonstrate function J = 5 N = 100 # Population item parameters as_t = rnorm(J,mean=2,sd=.5) bs_t = rnorm(J,mean=0,sd=.5) # Sampling gs and ss with truncation gs_t = rbeta(J,1,8) ps_g = pbeta(1-gs_t,1,8) ss_t = qbeta(runif(J)*ps_g,1,8) theta_t = rnorm(N) Y_t = Y_4pno_simulate(N,J,as=as_t,bs=bs_t,gs=gs_t,ss=ss_t,theta=theta_t) # Setting prior parameters mu_theta = 0 Sigma_theta_inv = 1 mu_xi = c(0,0) alpha_c = alpha_s = beta_c = beta_s = 1 Sigma_xi_inv = solve(2*matrix(c(1,0,0,1), 2, 2)) burnin = 1000 # Execute Gibbs sampler. This should take about 15.5 minutes out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,Sigma_theta_inv, alpha_c,beta_c,alpha_s, beta_s,burnin, rep(1,J),rep(1,J),gwg_reps=5,chain_length=burnin*2) # Summarizing posterior distribution OUT = cbind( apply(out_t$AS[, -c(1:burnin)], 1, mean), apply(out_t$BS[, -c(1:burnin)], 1, mean), apply(out_t$GS[, -c(1:burnin)], 1, mean), apply(out_t$SS[, -c(1:burnin)], 1, mean), apply(out_t$AS[, -c(1:burnin)], 1, sd), apply(out_t$BS[, -c(1:burnin)], 1, sd), apply(out_t$GS[, -c(1:burnin)], 1, sd), apply(out_t$SS[, -c(1:burnin)], 1, sd) ) OUT = cbind(1:J, OUT) colnames(OUT) = c('Item','as','bs','gs','ss','as_sd','bs_sd', 'gs_sd','ss_sd') print(OUT, digits = 3)#> Item as bs gs ss as_sd bs_sd gs_sd ss_sd #> [1,] 1 1.15 -0.107 0.230 0.2836 0.611 0.681 0.1317 0.1167 #> [2,] 2 1.11 0.353 0.284 0.2605 0.791 0.634 0.1173 0.1403 #> [3,] 3 1.39 -0.197 0.331 0.1859 0.658 0.861 0.1424 0.0990 #> [4,] 4 2.39 -0.457 0.240 0.0919 0.886 0.595 0.1403 0.0639 #> [5,] 5 2.30 0.167 0.133 0.1295 0.787 0.489 0.0757 0.0907