1 /******************************************************************
2 This file is a part of eco: R Package for Fitting Bayesian Models
3 of Ecological Inference for 2x2 Tables
4 by Kosuke Imai and Ying Lu
5 Copyright: GPL version 2 or later.
6 *******************************************************************/
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <Rmath.h>
12 #include <R_ext/Utils.h>
13 #include <R_ext/PrtUtil.h>
14 #include <R.h>
15 #include "vector.h"
16 #include "subroutines.h"
17 #include "rand.h"
18 #include "sample.h"
19 #include "macros.h"
20 #include "fintegrate.h"
21
22 /* Multivariate Normal density */
dMVN(double * Y,double * MEAN,double ** SIG_INV,int dim,int give_log)23 double dMVN(
24 double *Y, /* The data */
25 double *MEAN, /* The parameters */
26 double **SIG_INV, /* inverse of the covariance matrix */
27 int dim, /* dimension */
28 int give_log){ /* 1 if log_scale 0 otherwise */
29
30 int j,k;
31 double value=0.0;
32
33 for(j=0;j<dim;j++){
34 for(k=0;k<j;k++)
35 value+=2*(Y[k]-MEAN[k])*(Y[j]-MEAN[j])*SIG_INV[j][k];
36 value+=(Y[j]-MEAN[j])*(Y[j]-MEAN[j])*SIG_INV[j][j];
37 }
38
39 value=-0.5*value-0.5*dim*log(2*M_PI)+0.5*ddet(SIG_INV, dim, 1);
40
41
42 if(give_log)
43 return(value);
44 else
45 return(exp(value));
46
47 }
48
49
50 /* the density of Multivariate T-distribution */
dMVT(double * Y,double * MEAN,double ** SIG_INV,int nu,int dim,int give_log)51 double dMVT(
52 double *Y, /* The data */
53 double *MEAN, /* mean */
54 double **SIG_INV, /* inverse of scale matrix */
55 int nu, /* Degrees of freedom */
56 int dim, /* dimension */
57 int give_log) /* 1 if log_scale 0 otherwise */
58 {
59 int j,k;
60 double value=0;
61
62 for(j=0;j<dim;j++){
63 for(k=0;k<j;k++)
64 value+=2*(Y[k]-MEAN[k])*(Y[j]-MEAN[j])*SIG_INV[j][k];
65 value+=(Y[j]-MEAN[j])*(Y[j]-MEAN[j])*SIG_INV[j][j];
66 }
67
68 value=0.5*ddet(SIG_INV, dim,1) - 0.5*dim*(log((double)nu)+log(M_PI)) -
69 0.5*((double)dim+nu)*log(1+value/(double)nu) +
70 lgammafn(0.5*(double)(nu+dim)) - lgammafn(0.5*(double)nu);
71
72 if(give_log)
73 return(value);
74 else
75 return(exp(value));
76 }
77
78
79 /* Sample from the MVN dist */
rMVN(double * Sample,double * mean,double ** Var,int size)80 void rMVN(
81 double *Sample, /* Vector for the sample */
82 double *mean, /* The vector of means */
83 double **Var, /* The matrix Variance */
84 int size) /* The dimension */
85 {
86 int j,k;
87 double **Model = doubleMatrix(size+1, size+1);
88 double cond_mean;
89
90 /* draw from mult. normal using SWP */
91 for(j=1;j<=size;j++){
92 for(k=1;k<=size;k++)
93 Model[j][k]=Var[j-1][k-1];
94 Model[0][j]=mean[j-1];
95 Model[j][0]=mean[j-1];
96 }
97 Model[0][0]=-1;
98 Sample[0]=(double)norm_rand()*sqrt(Model[1][1])+Model[0][1];
99 for(j=2;j<=size;j++){
100 SWP(Model,j-1,size+1);
101 cond_mean=Model[j][0];
102 for(k=1;k<j;k++) cond_mean+=Sample[k-1]*Model[j][k];
103 Sample[j-1]=(double)norm_rand()*sqrt(Model[j][j])+cond_mean;
104 }
105
106 FreeMatrix(Model,size+1);
107 }
108
109
110 /* Sample from a wish dist */
111 /* Odell, P. L. and Feiveson, A. H. ``A Numerical Procedure to Generate
112 a Sample Covariance Matrix'' Journal of the American Statistical
113 Association, Vol. 61, No. 313. (Mar., 1966), pp. 199-203. */
114
rWish(double ** Sample,double ** S,int df,int size)115 void rWish(
116 double **Sample, /* The matrix with to hold the sample */
117 double **S, /* The parameter */
118 int df, /* the degrees of freedom */
119 int size) /* The dimension */
120 {
121 int i,j,k;
122 double *V = doubleArray(size);
123 double **B = doubleMatrix(size, size);
124 double **C = doubleMatrix(size, size);
125 double **N = doubleMatrix(size, size);
126 double **mtemp = doubleMatrix(size, size);
127
128 for(i=0;i<size;i++) {
129 V[i]=rchisq((double) df-i-1);
130 B[i][i]=V[i];
131 for(j=(i+1);j<size;j++)
132 N[i][j]=norm_rand();
133 }
134
135 for(i=0;i<size;i++) {
136 for(j=i;j<size;j++) {
137 Sample[i][j]=0;
138 Sample[j][i]=0;
139 mtemp[i][j]=0;
140 mtemp[j][i]=0;
141 if(i==j) {
142 if(i>0)
143 for(k=0;k<j;k++)
144 B[j][j]+=N[k][j]*N[k][j];
145 }
146 else {
147 B[i][j]=N[i][j]*sqrt(V[i]);
148 if(i>0)
149 for(k=0;k<i;k++)
150 B[i][j]+=N[k][i]*N[k][j];
151 }
152 B[j][i]=B[i][j];
153 }
154 }
155
156 dcholdc(S, size, C);
157 for(i=0;i<size;i++)
158 for(j=0;j<size;j++)
159 for(k=0;k<size;k++)
160 mtemp[i][j]+=C[i][k]*B[k][j];
161 for(i=0;i<size;i++)
162 for(j=0;j<size;j++)
163 for(k=0;k<size;k++)
164 Sample[i][j]+=mtemp[i][k]*C[j][k];
165
166 free(V);
167 FreeMatrix(B, size);
168 FreeMatrix(C, size);
169 FreeMatrix(N, size);
170 FreeMatrix(mtemp, size);
171 }
172
173 /* Sample from a Dirichlet distribution */
rDirich(double * Sample,double * theta,int size)174 void rDirich(
175 double *Sample, /* Vector for the sample */
176 double *theta, /* parameters */
177 int size) /* The dimension */
178 {
179 int j;
180 double dtemp=0;
181
182 for (j=0; j<size; j++) {
183 Sample[j] = rgamma(theta[j], 1.0);
184 dtemp += Sample[j];
185 }
186 for (j=0 ; j<size; j++)
187 Sample[j] /= dtemp;
188 }
189
190 /** density function on tomography line Y=XW_1+ (1-X)W_2
191 * Note: asssumes that the two points given W1* and W2*
192 * are on the tomography line
193 */
dBVNtomo(double * Wstar,void * pp,int give_log,double normc)194 double dBVNtomo(double *Wstar, /* Wstar values */
195 void* pp, //parameter
196 int give_log, /* 1 if log-scale, 0 otherwise */
197 double normc) //Normalization factor
198
199 {
200 int dim=2;
201 double *MEAN=doubleArray(dim);
202 double **SIGMA=doubleMatrix(dim,dim);
203 double density;
204 double rho, dtemp;
205
206 Param *param=(Param *)pp;
207 MEAN[0]=param->caseP.mu[0];
208 MEAN[1]=param->caseP.mu[1];
209 SIGMA[0][0]=param->setP->Sigma[0][0];
210 SIGMA[1][1]=param->setP->Sigma[1][1];
211 SIGMA[0][1]=param->setP->Sigma[0][1];
212 SIGMA[1][0]=param->setP->Sigma[1][0];
213
214
215 rho=SIGMA[0][1]/sqrt(SIGMA[0][0]*SIGMA[1][1]);
216 dtemp=1/(2*M_PI*sqrt(SIGMA[0][0]*SIGMA[1][1]*(1-rho*rho)));
217
218
219 density=-1/(2*(1-rho*rho))*
220 ((Wstar[0]-MEAN[0])*(Wstar[0]-MEAN[0])/SIGMA[0][0]+
221 +(Wstar[1]-MEAN[1])*(Wstar[1]-MEAN[1])/SIGMA[1][1]
222 -2*rho*(Wstar[0]-MEAN[0])*(Wstar[1]-MEAN[1])/sqrt(SIGMA[0][0]*SIGMA[1][1]))
223 +log(dtemp)-log(normc);
224
225 if (give_log==0) density=exp(density);
226 /*Rprintf("s11 %5g s22 %5g normc %5g dtemp %5g ldensity %5g\n", SIGMA[0][0],SIGMA[1][1],normc, dtemp, density);
227 char ch;
228 scanf(" %c", &ch );*/
229
230 Free(MEAN);
231 FreeMatrix(SIGMA,dim);
232
233 return density;
234
235
236 }
237
invLogit(double x)238 double invLogit(double x) {
239 if (x>30) return 0;
240 else return (1/(1+exp(-1*x)));
241 }
242
logit(double x,char * emsg)243 double logit(double x,char* emsg) {
244 if (x>=1 || x<=0) {
245 Rprintf(emsg);
246 Rprintf(": %5g is out of logit range\n",x);
247 }
248 return log(x/(1-x));
249 }
250
bit(int t,int n)251 int bit(int t, int n) {
252 t=t>>n;
253 return (t % 2);
254 }
255