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 /* Conditional Prediction for Normal Parametric Model for 2x2 Tables */
preBaseX(double * X,double * pdmu,double * pdSigma,int * pin_samp,int * pin_draw,int * verbose,double * pdStore)14 void preBaseX(
15 double *X, /* data X */
16 double *pdmu,
17 double *pdSigma,
18 int *pin_samp,
19 int *pin_draw,
20 int *verbose, /* 1 for output monitoring */
21 double *pdStore
22 ){
23
24 /* some integers */
25 int n_samp = *pin_samp; /* sample size */
26 int n_draw = *pin_draw; /* sample size of survey data */
27 int n_dim = 2;
28
29 double *mu = doubleArray(n_dim); /* The mean */
30 double *Wstar = doubleArray(n_dim);
31 double **Sigma = doubleMatrix(n_dim, n_dim); /* The covariance matrix */
32
33 /* misc variables */
34 int i, j, k, main_loop; /* used for various loops */
35 int itemp=0;
36 int itempM=0;
37 int itempS=0;
38 int progress = 1, itempP = ftrunc((double) n_draw/10);
39
40 /* get random seed */
41 GetRNGstate();
42
43 for(main_loop=0; main_loop<n_draw; main_loop++){
44 Sigma[0][0] = pdSigma[itempS]-pdSigma[itempS+2]*pdSigma[itempS+2]/pdSigma[itempS+5];
45 Sigma[1][1] = pdSigma[itempS+3]-pdSigma[itempS+4]*pdSigma[itempS+4]/pdSigma[itempS+5];
46 Sigma[0][1] = pdSigma[itempS+1]-pdSigma[itempS+2]*pdSigma[itempS+4]/pdSigma[itempS+5];
47 Sigma[1][0] = Sigma[0][1];
48 for(i=0; i<n_samp; i++) {
49 mu[0] = pdmu[itempM]+pdSigma[itempS+2]/pdSigma[itempS+5]*(X[i]-pdmu[itempM+2]);
50 mu[1] = pdmu[itempM+1]+pdSigma[itempS+4]/pdSigma[itempS+5]*(X[i]-pdmu[itempM+2]);
51 rMVN(Wstar, mu, Sigma, n_dim);
52 for (j=0; j<n_dim; j++)
53 pdStore[itemp++] = exp(Wstar[j])/(1+exp(Wstar[j]));
54 }
55 itempS += 6;
56 itempM += 3;
57 if (*verbose)
58 if (itempP == main_loop) {
59 Rprintf("%3d percent done.\n", progress*10);
60 itempP+=ftrunc((double) n_draw/10); progress++;
61 R_FlushConsole();
62 }
63 R_CheckUserInterrupt();
64 }
65
66 if(*verbose)
67 Rprintf("100 percent done.\n");
68
69 /** write out the random seed **/
70 PutRNGstate();
71
72 /* Freeing the memory */
73 free(mu);
74 free(Wstar);
75 FreeMatrix(Sigma,n_dim);
76
77 } /* main */
78
79