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.h>
14 #include "vector.h"
15 #include "subroutines.h"
16 #include "rand.h"
17
18 /** Normal-InvWishart updating
19 Y|mu, Sigma ~ N(mu, Sigma)
20 mu|Sigma ~ N(mu0, Sigma/tau0)
21 Sigma ~ InvWish(nu0, S0^{-1}) **/
NIWupdate(double ** Y,double * mu,double ** Sigma,double ** InvSigma,double * mu0,double tau0,int nu0,double ** S0,int n_samp,int n_dim)22 void NIWupdate(
23 double **Y, /* data */
24 double *mu, /* mean */
25 double **Sigma, /* variance */
26 double **InvSigma, /* precision */
27 double *mu0, /* prior mean */
28 double tau0, /* prior scale */
29 int nu0, /* prior df */
30 double **S0, /* prior scale */
31 int n_samp, /* sample size */
32 int n_dim) /* dimension */
33 {
34 int i,j,k;
35 double *Ybar = doubleArray(n_dim);
36 double *mun = doubleArray(n_dim);
37 double **Sn = doubleMatrix(n_dim, n_dim);
38 double **mtemp = doubleMatrix(n_dim, n_dim);
39
40 /*read data */
41 for (j=0; j<n_dim; j++) {
42 Ybar[j] = 0;
43 for (i=0; i<n_samp; i++)
44 Ybar[j] += Y[i][j];
45 Ybar[j] /= n_samp;
46 for (k=0; k<n_dim; k++)
47 Sn[j][k] = S0[j][k];
48 }
49
50 /* posterior updating*/
51
52 for (j=0; j<n_dim; j++)
53 {
54 mun[j] = (tau0*mu0[j]+n_samp*Ybar[j])/(tau0+n_samp);
55 for (k=0; k<n_dim; k++)
56 {
57 Sn[j][k] += (tau0*n_samp)*(Ybar[j]-mu0[j])*(Ybar[k]-mu0[k])/(tau0+n_samp);
58 for (i=0; i<n_samp; i++)
59 Sn[j][k] += (Y[i][j]-Ybar[j])*(Y[i][k]-Ybar[k]);
60 }
61 }
62
63 dinv(Sn, n_dim, mtemp);
64 rWish(InvSigma, mtemp, nu0+n_samp, n_dim);
65 dinv(InvSigma, n_dim, Sigma);
66
67 for (j=0; j<n_dim; j++)
68 for (k=0; k<n_dim; k++)
69 mtemp[j][k] = Sigma[j][k]/(tau0+n_samp);
70
71 rMVN(mu, mun, mtemp, n_dim);
72
73 free(Ybar);
74 free(mun);
75 FreeMatrix(Sn, n_dim);
76 FreeMatrix(mtemp, n_dim);
77 }
78