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