1 #include <string.h>
2 #include <stddef.h>
3 #include <stdio.h>
4 #include <math.h>
5 #include <Rmath.h>
6 #include <R.h>
7 #include "vector.h"
8 #include "subroutines.h"
9 #include "rand.h"
10 #include "bayes.h"
11 #include "sample.h"
12
13 /* Normal Parametric Model for 2x2 Tables with Contextual Effects */
cBaseecoX(double * pdX,int * pin_samp,int * n_gen,int * burn_in,int * pinth,int * verbose,int * pinu0,double * pdtau0,double * mu0,double * pdS0,double * mustart,double * Sigmastart,int * survey,int * sur_samp,double * sur_W,int * x1,int * sampx1,double * x1_W1,int * x0,int * sampx0,double * x0_W2,double * minW1,double * maxW1,int * parameter,int * Grid,double * pdSMu0,double * pdSMu1,double * pdSMu2,double * pdSSig00,double * pdSSig01,double * pdSSig02,double * pdSSig11,double * pdSSig12,double * pdSSig22,double * pdSW1,double * pdSW2)14 void cBaseecoX(
15 /*data input */
16 double *pdX, /* data (X, Y) */
17 int *pin_samp, /* sample size */
18
19 /*MCMC draws */
20 int *n_gen, /* number of gibbs draws */
21 int *burn_in, /* number of draws to be burned in */
22 int *pinth, /* keep every nth draw */
23 int *verbose, /* 1 for output monitoring */
24
25 /* prior specification*/
26 int *pinu0, /* prior df parameter for InvWish */
27 double *pdtau0, /* prior scale parameter for Sigma under G0*/
28 double *mu0, /* prior mean for mu under G0 */
29 double *pdS0, /* prior scale for Sigma */
30 double *mustart, /* starting values for mu */
31 double *Sigmastart, /* starting values for Sigma */
32
33 /*incorporating survey data */
34 int *survey, /*1 if survey data available (set of W_1, W_2)
35 0 not*/
36 int *sur_samp, /*sample size of survey data*/
37 double *sur_W, /*set of known W_1, W_2 */
38
39 /* incorporating homeogenous areas */
40 int *x1, /* 1 if X=1 type areas available W_1 known, W_2 unknown */
41 int *sampx1, /* number X=1 type areas */
42 double *x1_W1, /* values of W_1 for X1 type areas */
43
44 int *x0, /* 1 if X=0 type areas available W_2
45 known, W_1 unknown */
46 int *sampx0, /* number X=0 type areas */
47 double *x0_W2, /* values of W_2 for X0 type areas */
48
49 /* bounds fo W1 */
50 double *minW1, double *maxW1,
51
52 /* flags */
53 int *parameter, /* 1 if save population parameter */
54 int *Grid, /* 1 if Grid algorithm is used; 0 for
55 Metropolis */
56
57 /* storage for Gibbs draws of mu/sigmat*/
58 double *pdSMu0, double *pdSMu1, double *pdSMu2,
59 double *pdSSig00, double *pdSSig01, double *pdSSig02,
60 double *pdSSig11, double *pdSSig12, double *pdSSig22,
61
62 /* storage for Gibbs draws of W*/
63 double *pdSW1, double *pdSW2
64 ){
65
66 /* some integers */
67 int n_samp = *pin_samp; /* sample size */
68 int s_samp = *sur_samp; /* sample size of survey data */
69 int x1_samp = *sampx1; /* sample size for X=1 */
70 int x0_samp = *sampx0; /* sample size for X=0 */
71 int t_samp = n_samp+s_samp+x1_samp+x0_samp; /* total sample size */
72 int nth = *pinth;
73 int n_dim = 2; /* dimension */
74 int n_step = 1000; /* 1/The default size of grid step */
75
76 /* prior parameters */
77 double tau0 = *pdtau0;
78 int nu0 = *pinu0;
79 double **S0 = doubleMatrix(n_dim+1,n_dim+1); /* The prior S parameter for InvWish */
80
81 /* data */
82 double **X = doubleMatrix(n_samp,n_dim); /* The Y and covariates */
83 double **W = doubleMatrix(t_samp,n_dim); /* The W1 and W2 matrix */
84 double **Wstar = doubleMatrix(t_samp,n_dim+1); /* logit transformed
85 W and X */
86 double **S_W = doubleMatrix(s_samp, n_dim+1); /* known W1, W2, X */
87 double **S_Wstar = doubleMatrix(s_samp, n_dim+1); /* logit
88 transformed S_W */
89 /* grids */
90 double **W1g = doubleMatrix(n_samp, n_step);
91 double **W2g = doubleMatrix(n_samp, n_step);
92 int *n_grid = intArray(n_samp); /* grid size */
93
94 /* ordinary model variables */
95 double *mu = doubleArray(n_dim+1);
96 double **Sigma = doubleMatrix(n_dim+1,n_dim+1);
97 double **InvSigma = doubleMatrix(n_dim+1,n_dim+1);
98
99 /* conditional mean & variance for (W1, W2) given X */
100 double *mu_w = doubleArray(n_dim);
101 double **Sigma_w = doubleMatrix(n_dim,n_dim);
102 double **InvSigma_w = doubleMatrix(n_dim,n_dim);
103
104 /* misc variables */
105 int i, j, k, t, main_loop; /* used for various loops */
106 int itemp, itempS, itempC, itempA;
107 int progress = 1, itempP = ftrunc((double) *n_gen/10);
108 double dtemp, dtemp1;
109
110 /* get random seed */
111 GetRNGstate();
112
113 /* priors */
114 itemp = 0;
115 for(k=0; k<(n_dim+1); k++)
116 for(j=0; j<(n_dim+1); j++)
117 S0[j][k] = pdS0[itemp++];
118
119 /* read the data set */
120 itemp = 0;
121 for (j = 0; j < n_dim; j++)
122 for (i = 0; i < n_samp; i++)
123 X[i][j] = pdX[itemp++];
124
125 /* Initialize W, Wstar for n_samp */
126 for (i=0; i< n_samp; i++) {
127 if (X[i][1]!=0 && X[i][1]!=1) {
128 W[i][0]=runif(minW1[i], maxW1[i]);
129 W[i][1]=(X[i][1]-X[i][0]*W[i][0])/(1-X[i][0]);
130 }
131 if (X[i][1]==0)
132 for (j=0; j<n_dim; j++) W[i][j]=0.0001;
133 if (X[i][1]==1)
134 for (j=0; j<n_dim; j++) W[i][j]=0.9999;
135 for (j=0; j<n_dim; j++)
136 Wstar[i][j]=log(W[i][j])-log(1-W[i][j]);
137 Wstar[i][n_dim] = log(X[i][0])-log(1-X[i][0]);
138 }
139
140 /*read homeogenous areas information */
141 if (*x1==1)
142 for (i=0; i<x1_samp; i++) {
143 W[(n_samp+i)][0]=x1_W1[i];
144 if (W[(n_samp+i)][0]==0) W[(n_samp+i)][0]=0.0001;
145 if (W[(n_samp+i)][0]==1) W[(n_samp+i)][0]=0.9999;
146 Wstar[(n_samp+i)][0]=log(W[(n_samp+i)][0])-log(1-W[(n_samp+i)][0]);
147 Wstar[n_samp+i][n_dim]=log(0.9999)-log(0.0001);
148 }
149
150 if (*x0==1)
151 for (i=0; i<x0_samp; i++) {
152 W[(n_samp+x1_samp+i)][1]=x0_W2[i];
153 if (W[(n_samp+x1_samp+i)][1]==0) W[(n_samp+x1_samp+i)][1]=0.0001;
154 if (W[(n_samp+x1_samp+i)][1]==1) W[(n_samp+x1_samp+i)][1]=0.9999;
155 Wstar[(n_samp+x1_samp+i)][1]=log(W[(n_samp+x1_samp+i)][1])-log(1-W[(n_samp+x1_samp+i)][1]);
156 Wstar[n_samp+x1_samp+i][n_dim]=log(0.0001)-log(0.9999);
157 }
158
159 /*read the survey data */
160 if (*survey==1) {
161 itemp = 0;
162 for (j=0; j<=n_dim; j++)
163 for (i=0; i<s_samp; i++) {
164 S_W[i][j]=sur_W[itemp++];
165 if (S_W[i][j]==0) S_W[i][j]=0.0001;
166 if (S_W[i][j]==1) S_W[i][j]=0.9999;
167 S_Wstar[i][j]=log(S_W[i][j])-log(1-S_W[i][j]);
168 if (j<n_dim) {
169 W[(n_samp+x1_samp+x0_samp+i)][j]=S_W[i][j];
170 Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
171 }
172 else
173 Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
174 }
175 }
176
177 /* counters */
178 itempA=0; /* for alpha */
179 itempS=0; /* for storage */
180 itempC=0; /* control nth draw */
181
182 /*** calculate grids ***/
183 if (*Grid)
184 GridPrep(W1g, W2g, X, maxW1, minW1, n_grid, n_samp, n_step);
185
186 /* starting values of mu and Sigma */
187 itemp = 0;
188 for(j=0;j<(n_dim+1);j++){
189 mu[j] = mustart[j];
190 for(k=0;k<(n_dim+1);k++)
191 Sigma[j][k]=Sigmastart[itemp++];
192 }
193 dinv(Sigma, n_dim+1, InvSigma);
194
195 /***Gibbs Sampler ***/
196 if (*verbose)
197 Rprintf("Starting Gibbs Sampler...\n");
198 for(main_loop=0; main_loop<*n_gen; main_loop++){
199 /* conditional variance */
200 for (j=0; j<n_dim; j++)
201 for (k=0; k<n_dim; k++)
202 Sigma_w[j][k]=Sigma[j][k]-Sigma[n_dim][j]/Sigma[n_dim][n_dim]*Sigma[n_dim][k];
203 dinv(Sigma_w, n_dim, InvSigma_w);
204
205 /**update W, Wstar given mu, Sigma in regular areas**/
206 for (i=0; i<n_samp; i++){
207 for (j=0; j<n_dim; j++)
208 mu_w[j]=mu[j]+Sigma[n_dim][j]/Sigma[n_dim][n_dim]*(Wstar[i][2]-mu[n_dim]);
209 if ( X[i][1]!=0 && X[i][1]!=1 ) {
210 if (*Grid)
211 rGrid(W[i], W1g[i],W2g[i], n_grid[i], mu_w, InvSigma_w,
212 n_dim);
213 else
214 rMH(W[i], X[i], minW1[i], maxW1[i], mu_w, InvSigma_w, n_dim);
215 }
216 /*3 compute Wsta_i from W_i*/
217 Wstar[i][0]=log(W[i][0])-log(1-W[i][0]);
218 Wstar[i][1]=log(W[i][1])-log(1-W[i][1]);
219 }
220
221 /*update W2 given W1, mu and Sigma in x1 homeogeneous areas */
222 if (*x1==1)
223 for (i=0; i<x1_samp; i++) {
224 dtemp=mu_w[1]+Sigma_w[0][1]/Sigma_w[0][0]*(Wstar[n_samp+i][0]-mu_w[0]);
225 dtemp1=Sigma_w[1][1]*(1-Sigma_w[0][1]*Sigma_w[0][1]/(Sigma_w[0][0]*Sigma_w[1][1]));
226 dtemp1=sqrt(dtemp1);
227 Wstar[n_samp+i][1]=rnorm(dtemp, dtemp1);
228 W[n_samp+i][1]=exp(Wstar[n_samp+i][1])/(1+exp(Wstar[n_samp+i][1]));
229 }
230
231 /*update W1 given W2, mu and Sigma in x0 homeogeneous areas */
232 if (*x0==1)
233 for (i=0; i<x0_samp; i++) {
234 dtemp=mu_w[0]+Sigma_w[0][1]/Sigma_w[1][1]*(Wstar[n_samp+x1_samp+i][1]-mu_w[1]);
235 dtemp1=Sigma_w[0][0]*(1-Sigma_w[0][1]*Sigma_w[0][1]/(Sigma_w[0][0]*Sigma_w[1][1]));
236 dtemp1=sqrt(dtemp1);
237 Wstar[n_samp+x1_samp+i][0]=rnorm(dtemp, dtemp1);
238 W[n_samp+x1_samp+i][0]=exp(Wstar[n_samp+x1_samp+i][0])/(1+exp(Wstar[n_samp+x1_samp+i][0]));
239 }
240
241 /* update mu, Sigma given wstar using effective sample of Wstar */
242 NIWupdate(Wstar, mu, Sigma, InvSigma, mu0, tau0, nu0, S0, t_samp, n_dim+1);
243
244 /*store Gibbs draw after burn-in and every nth draws */
245 R_CheckUserInterrupt();
246 if (main_loop>=*burn_in){
247 itempC++;
248 if (itempC==nth){
249 pdSMu0[itempA]=mu[0];
250 pdSMu1[itempA]=mu[1];
251 pdSMu2[itempA]=mu[2];
252 pdSSig00[itempA]=Sigma[0][0];
253 pdSSig01[itempA]=Sigma[0][1];
254 pdSSig02[itempA]=Sigma[0][2];
255 pdSSig11[itempA]=Sigma[1][1];
256 pdSSig12[itempA]=Sigma[1][2];
257 pdSSig22[itempA]=Sigma[2][2];
258 itempA++;
259 for(i=0; i<(n_samp+x1_samp+x0_samp); i++){
260 pdSW1[itempS]=W[i][0];
261 pdSW2[itempS]=W[i][1];
262 itempS++;
263 }
264 itempC=0;
265 }
266 } /*end of stroage *burn_in*/
267 if (*verbose)
268 if (itempP == main_loop) {
269 Rprintf("%3d percent done.\n", progress*10);
270 itempP+=ftrunc((double) *n_gen/10); progress++;
271 R_FlushConsole();
272 }
273 } /*end of MCMC for normal */
274
275 if(*verbose)
276 Rprintf("100 percent done.\n");
277
278
279 /** write out the random seed **/
280 PutRNGstate();
281
282 /* Freeing the memory */
283 FreeMatrix(X, n_samp);
284 FreeMatrix(W, t_samp);
285 FreeMatrix(Wstar, t_samp);
286 free(n_grid);
287 FreeMatrix(S0, n_dim+1);
288 FreeMatrix(W1g, n_samp);
289 FreeMatrix(W2g, n_samp);
290 FreeMatrix(S_W, s_samp);
291 FreeMatrix(S_Wstar, s_samp);
292 free(mu);
293 FreeMatrix(Sigma, n_dim+1);
294 FreeMatrix(InvSigma, n_dim+1);
295 free(mu_w);
296 FreeMatrix(Sigma_w, n_dim);
297 FreeMatrix(InvSigma_w, n_dim);
298 } /* main */
299
300