1 #include <stddef.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <Rmath.h>
5 #include <R.h>
6 #include "vector.h"
7 #include "subroutines.h"
8 #include "rand.h"
9 #include "bayes.h"
10 #include "sample.h"
11
12 /* Normal Parametric Model for 2xC (with C > 2) Tables */
cBase2C(double * pdX,double * Y,double * pdWmin,double * pdWmax,int * pin_samp,int * pin_col,int * reject,int * maxit,int * n_gen,int * burn_in,int * pinth,int * verbose,int * pinu0,double * pdtau0,double * mu0,double * pdS0,double * mu,double * SigmaStart,int * parameter,double * pdSmu,double * pdSSigma,double * pdSW)13 void cBase2C(
14 /*data input */
15 double *pdX, /* X: matrix */
16 double *Y, /* Y: vector */
17 double *pdWmin, /* lower bounds */
18 double *pdWmax, /* uppwer bounds */
19 int *pin_samp, /* sample size */
20 int *pin_col, /* number of columns */
21
22 /*MCMC draws */
23 int *reject, /* whether to use rejection sampling */
24 int *maxit, /* max number of iterations for
25 rejection sampling */
26 int *n_gen, /* number of gibbs draws */
27 int *burn_in, /* number of draws to be burned in */
28 int *pinth, /* keep every nth draw */
29 int *verbose, /* 1 for output monitoring */
30
31 /* prior specification*/
32 int *pinu0, /* prior df parameter for InvWish */
33 double *pdtau0, /* prior scale parameter for Sigma */
34 double *mu0, /* prior mean for mu */
35 double *pdS0, /* prior scale for Sigma */
36
37 /* starting values */
38 double *mu,
39 double *SigmaStart,
40
41 /* storage */
42 int *parameter, /* 1 if save population parameter */
43 double *pdSmu,
44 double *pdSSigma,
45 double *pdSW
46 ){
47
48 /* some integers */
49 int n_samp = *pin_samp; /* sample size */
50 int nth = *pinth; /* keep every pth draw */
51 int n_col = *pin_col; /* dimension */
52
53 /* prior parameters */
54 double tau0 = *pdtau0; /* prior scale for mu */
55 int nu0 = *pinu0; /* prior degrees of freedom */
56 double **S0 = doubleMatrix(n_col, n_col); /* prior scale for Sigma */
57
58 /* data */
59 double **X = doubleMatrix(n_samp, n_col); /* X */
60 double **W = doubleMatrix(n_samp, n_col); /* The W matrix */
61 double **Wstar = doubleMatrix(n_samp, n_col); /* logit(W) */
62
63 /* The lower and upper bounds of U = W*X/Y **/
64 double **minU = doubleMatrix(n_samp, n_col);
65 double **maxU = doubleMatrix(n_samp, n_col);
66
67 /* model parameters */
68 double **Sigma = doubleMatrix(n_col, n_col); /* The covariance matrix */
69 double **InvSigma = doubleMatrix(n_col, n_col); /* The inverse covariance matrix */
70
71 /* misc variables */
72 int i, j, k, main_loop; /* used for various loops */
73 int itemp;
74 int itempM = 0; /* for mu */
75 int itempS = 0; /* for Sigma */
76 int itempW = 0; /* for W */
77 int itempC = 0; /* control nth draw */
78 int progress = 1, itempP = ftrunc((double) *n_gen/10);
79 double dtemp, dtemp1;
80 double *param = doubleArray(n_col); /* Dirichlet parameters */
81 double *dvtemp = doubleArray(n_col);
82
83 /* get random seed */
84 GetRNGstate();
85
86 /* read X */
87 itemp = 0;
88 for (j = 0; j < n_col; j++)
89 for (i = 0; i < n_samp; i++)
90 X[i][j] = pdX[itemp++];
91
92 /* read initial values of Sigma */
93 itemp = 0;
94 for (k = 0; k < n_col; k++)
95 for (j = 0; j < n_col; j++)
96 Sigma[j][k] = SigmaStart[itemp++];
97 dinv(Sigma, n_col, InvSigma);
98
99 /* compute bounds on U */
100 itemp = 0;
101 for (j = 0; j < n_col; j++)
102 for (i = 0; i < n_samp; i++)
103 minU[i][j] = fmax2(0, pdWmin[itemp++]*X[i][j]/Y[i]);
104 itemp = 0;
105 for (j = 0; j < n_col; j++)
106 for (i = 0; i < n_samp; i++)
107 maxU[i][j] = fmin2(1, pdWmax[itemp++]*X[i][j]/Y[i]);
108
109 /* initial values for W */
110 for (j = 0; j < n_col; j++)
111 param[j] = 1;
112 for (i = 0; i < n_samp; i++) {
113 k = 0; itemp = 1;
114 while (itemp > 0) { /* rejection sampling */
115 rDirich(dvtemp, param, n_col);
116 itemp = 0; k++;
117 for (j = 0; j < n_col; j++)
118 if (dvtemp[j] > maxU[i][j] || dvtemp[j] < minU[i][j])
119 itemp++;
120 if (itemp == 0)
121 for (j = 0; j < n_col; j++)
122 W[i][j] = dvtemp[j]*Y[i]/X[i][j];
123 if (k > *maxit) { /* if rejection sampling fails, then use
124 midpoits of bounds sequentially */
125 itemp = 0;
126 dtemp = Y[i]; dtemp1 = 1;
127 for (j = 0; j < n_col-1; j++) {
128 W[i][j] = 0.5*(fmax2(0,(X[i][j]/dtemp1+dtemp-1)*dtemp1/X[i][j])+
129 fmin2(1,dtemp*dtemp1/X[i][j]));
130 dtemp -= W[i][j]*X[i][j]/dtemp1;
131 dtemp1 -= X[i][j];
132 }
133 W[i][n_col-1] = dtemp;
134 }
135 }
136 for (j = 0; j < n_col; j++)
137 Wstar[i][j] = log(W[i][j])-log(1-W[i][j]);
138 }
139
140 /* read the prior */
141 itemp = 0;
142 for(k = 0; k < n_col; k++)
143 for(j = 0; j < n_col; j++)
144 S0[j][k] = pdS0[itemp++];
145
146 /*** Gibbs sampler! ***/
147 if (*verbose)
148 Rprintf("Starting Gibbs sampler...\n");
149 for(main_loop = 0; main_loop < *n_gen; main_loop++){
150 /** update W, Wstar given mu, Sigma **/
151 for (i = 0; i < n_samp; i++){
152 rMH2c(W[i], X[i], Y[i], minU[i], maxU[i], mu, InvSigma, n_col,
153 *maxit, *reject);
154 for (j = 0; j < n_col; j++)
155 Wstar[i][j] = log(W[i][j])-log(1-W[i][j]);
156 }
157
158 /* update mu, Sigma given wstar using effective sample of Wstar */
159 NIWupdate(Wstar, mu, Sigma, InvSigma, mu0, tau0, nu0, S0, n_samp, n_col);
160
161 /*store Gibbs draw after burn-in and every nth draws */
162 if (main_loop>=*burn_in){
163 itempC++;
164 if (itempC==nth){
165 for (j = 0; j < n_col; j++) {
166 pdSmu[itempM++]=mu[j];
167 for (k = 0; k < n_col; k++)
168 if (j <=k)
169 pdSSigma[itempS++]=Sigma[j][k];
170 }
171 for(i = 0; i < n_samp; i++)
172 for (j = 0; j < n_col; j++)
173 pdSW[itempW++] = W[i][j];
174 itempC=0;
175 }
176 }
177 if (*verbose)
178 if (itempP == main_loop) {
179 Rprintf("%3d percent done.\n", progress*10);
180 itempP+=ftrunc((double) *n_gen/10); progress++;
181 R_FlushConsole();
182 }
183 R_CheckUserInterrupt();
184 } /* end of Gibbs sampler */
185
186 if(*verbose)
187 Rprintf("100 percent done.\n");
188
189 /** write out the random seed **/
190 PutRNGstate();
191
192 /* Freeing the memory */
193 FreeMatrix(S0, n_col);
194 FreeMatrix(X, n_samp);
195 FreeMatrix(W, n_samp);
196 FreeMatrix(Wstar, n_samp);
197 FreeMatrix(minU, n_samp);
198 FreeMatrix(maxU, n_samp);
199 FreeMatrix(Sigma, n_col);
200 FreeMatrix(InvSigma, n_col);
201 free(dvtemp);
202 free(param);
203 } /* main */
204
205