diff options
Diffstat (limited to 'sourcecodes/localscore/modified_postc.c')
| -rw-r--r-- | sourcecodes/localscore/modified_postc.c | 171 |
1 files changed, 171 insertions, 0 deletions
diff --git a/sourcecodes/localscore/modified_postc.c b/sourcecodes/localscore/modified_postc.c new file mode 100644 index 00000000..f223fdb8 --- /dev/null +++ b/sourcecodes/localscore/modified_postc.c @@ -0,0 +1,171 @@ +/* -*- Mode: C -*- + * postc.c --- Posterior for continuous node with continuous parents + * Author : Claus Dethlefsen + * Created On : Tue Mar 12 06:44:35 2002 + * Last Modified By: Claus Dethlefsen + * Last Modified On: Wed Jun 04 11:56:51 2003 + * Update Count : 227 + * Status : Unknown, Use with caution! + */ + +/* + ## +## Copyright (C) 2002 Susanne Gammelgaard Bøttcher, Claus Dethlefsen +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +###################################################################### +*/ + +#include "R.h" +#include "Rmath.h" +#include "matrix.h" +#include "modified_matrix.c" + +void postc(double *mu, double *tau, double *rho, double *phi, double + *loglik, double *y, double *z, int *n, int *d) +{ + int i,j,ii,jj; + double logscale,logk,mscore; + double **mtau, **mmu, **tauinv=0; + double **zero, **zi, **ziy; + double **oldtau, **oldmu; + //temp pointers + double **tm1,**tm2,**tm3,**tm4,**mm1,**sm1,**sm2,**tt1,**tt2; + + /* allocate space for matrices */ + mtau = dmatrix(1,*d,1,*d); + oldtau = dmatrix(1,*d,1,*d); + + zi = dmatrix(1,*d,1,1); + ziy = dmatrix(1,*d,1,1); + mmu = dmatrix(1,*d,1,1); + oldmu = dmatrix(1,*d,1,1); + zero = dmatrix(1,*d,1,1); + tauinv = dmatrix(1,*d,1,*d); + + + + /* copy arguments into the matrices */ + asmatrix(mu,mmu,*d,1); + asmatrix(tau,mtau,*d,*d); + + /* show input */ + + for(i = 1; i <= *n; i++) { + for (ii=1; ii<=*d; ii++) { + for (jj=1; jj<=*d; jj++) { + tauinv[ii][jj] = mtau[ii][jj]; + } + } + + invers(tauinv, *d, zero, 1); + + for (j=1; j<=*d; j++) { + zi[j][1] = z[j-1+(i-1)*(*d)]; + } + //define once + tt1=transp(zi,*d,1); + + tm1=matmult(tauinv,zi,*d,*d,1); + tm2=matmult(tt1,tm1,1,*d,1); + + logscale = log(*phi) + log1p(tm2[1][1]); + + free_dmatrix(tm1,1,*d,1,*d); + free_dmatrix(tm2,1,1,1,*d); + + + logk = lgammafn( 0.5*(1.0+*rho) ) - lgammafn(*rho*0.5); + logk -= 0.5*(logscale + log(M_PI)); + + tm1=matmult(tt1,mmu,1,*d,1); + + + mscore = logk - 0.5*(*rho+1)*log1p((y[i-1] - tm1[1][1])*(y[i-1] - tm1[1][1])/exp(logscale)); + *loglik += mscore; + + free_dmatrix(tm1,1,1,1,1); + + for (ii=1; ii<=*d; ii++) { + for (jj=1; jj<=*d; jj++) { + oldtau[ii][jj] = mtau[ii][jj]; + } + } + + + for (jj=1; jj<=*d; jj++) { + oldmu[jj][1] = mmu[jj][1]; + } + + tm1=matmult(zi,tt1,*d,1,*d); + sm1=matsum(mtau,tm1,*d,*d); + free_dmatrix(mtau,1,*d,1,*d); + mtau = sm1; + free_dmatrix(tm1,1,*d,1,*d); + + for (ii=1; ii<=*d; ii++) { + for (jj=1; jj<=*d; jj++) { + tauinv[ii][jj] = mtau[ii][jj]; + } + } + + + invers(tauinv, *d, zero, 1); + + for (j=1;j<=*d;j++) + ziy[j][1] = zi[j][1]*y[i-1]; + + + tm1=matmult(oldtau,mmu,*d,*d,1); + sm2=matsum(tm1,ziy,*d,1); + free_dmatrix(mmu,1,*d,1,1); + mmu=matmult(tauinv,sm2,*d,*d,1); + free_dmatrix(tm1,1,*d,1,*d); + free_dmatrix(sm2,1,*d,1,1); + + + (*rho)++; + + tm1=matmult(tt1,mmu,1,*d,1); + mm1=matminus(oldmu,mmu,*d,1); + tt2=transp(mm1,*d,1); + tm3=matmult(oldtau,oldmu,*d,*d,1); + tm4=matmult(tt2,tm3,1,*d,1); + + (*phi) += (y[i-1]-tm1[1][1])*y[i-1] + tm4[1][1]; + + + free_dmatrix(tm1,1,1,1,*d); + free_dmatrix(mm1,1,*d,1,1); + free_dmatrix(tt2,1,1,1,*d); + free_dmatrix(tm3,1,*d,1,*d); + free_dmatrix(tm4,1,1,1,*d); + + free_dmatrix(tt1,1,1,1,*d); + + + } + + free_dmatrix(mtau,1,*d,1,*d); + free_dmatrix(zi,1,*d,1,1); + free_dmatrix(ziy,1,*d,1,1); + free_dmatrix(mmu,1,*d,1,1); + free_dmatrix(zero,1,*d,1,1); + free_dmatrix(tauinv,1,*d,1,*d); + free_dmatrix(oldtau,1,*d,1,*d); + free_dmatrix(oldmu,1,*d,1,1); + +} + |