[Depmix-commits] r468 - in pkg/depmixS4: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 20 14:42:28 CEST 2011
Author: ingmarvisser
Date: 2011-06-20 14:42:27 +0200 (Mon, 20 Jun 2011)
New Revision: 468
Modified:
pkg/depmixS4/R/fb.R
pkg/depmixS4/src/fb.cc
pkg/depmixS4/src/fb.h
Log:
Loop over cases and times for forward variables now complete (still to add for the mix case, ie no transition model).
Modified: pkg/depmixS4/R/fb.R
===================================================================
--- pkg/depmixS4/R/fb.R 2011-06-20 08:01:49 UTC (rev 467)
+++ pkg/depmixS4/R/fb.R 2011-06-20 12:42:27 UTC (rev 468)
@@ -24,34 +24,47 @@
nt <- nrow(B)
ns <- ncol(init)
- alpha <- matrix(0,ncol=ns,nrow=nt)
- beta <- matrix(0,ncol=ns,nrow=nt)
- sca <- rep(0,nt)
- xi <- array(0,dim=c(nt,ns,ns))
-
if(is.null(ntimes)) ntimes <- nt
lt <- length(ntimes)
-
- if(useC) {
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ if(useC) {
+
+ alpha <- matrix(0,ncol=ns,nrow=nt)
+ beta <- matrix(0,ncol=ns,nrow=nt)
+ sca <- rep(0,nt)
+ xi <- array(0,dim=c(nt,ns,ns))
+
res <- .C("forwardbackward",
- as.integer(ns),
- as.integer(lt),
- as.integer(nt),
- as.integer(ntimes),
- as.double(init),
- as.double(A),
- as.double(B),
- as.double(alpha),
- as.double(beta),
- as.double(sca),
- as.double(xi),
- package="depmixS4")
+ ns=as.integer(ns),
+ lt=as.integer(lt),
+ nt=as.integer(nt),
+ ntimes=as.integer(ntimes),
+ bt=as.integer(bt),
+ et=as.integer(et),
+ init=as.double(t(init)),
+ A=as.double(A),
+ B=as.double(t(B)),
+ alpha=as.double(alpha),
+ beta=as.double(beta),
+ sca=as.double(sca),
+ xi=as.double(xi),
+ package="depmixS4")[c("alpha","beta","sca","xi")]
+ res$alpha <- matrix(res$alpha,nc=ns,byrow=TRUE)
+ res$beta <- matrix(res$beta,nc=ns,byrow=TRUE)
+
+ res$gamma <- alpha*beta/sca
+ res$logLike <- -sum(log(sca))
+
} else {
- et <- cumsum(ntimes)
- bt <- c(1,et[-lt]+1)
+ alpha <- matrix(ncol=ns,nrow=nt)
+ beta <- matrix(ncol=ns,nrow=nt)
+ sca <- vector(length=nt)
+ xi <- array(dim=c(nt,ns,ns))
for(case in 1:lt) {
alpha[bt[case],] <- init[case,]*B[bt[case],] # initialize
@@ -91,7 +104,7 @@
} else {
res <- list(gamma=gamma,xi=xi,logLike=like)
}
- }
+ }
res
}
Modified: pkg/depmixS4/src/fb.cc
===================================================================
--- pkg/depmixS4/src/fb.cc 2011-06-20 08:01:49 UTC (rev 467)
+++ pkg/depmixS4/src/fb.cc 2011-06-20 12:42:27 UTC (rev 468)
@@ -41,7 +41,7 @@
// gamma is computed as alpha*beta/sca in R (no loop needed)
-void forwardbackward(int *ns, int *nc, int *nt, int *ntimes,
+void forwardbackward(int *ns, int *nc, int *nt, int *ntimes, int *bt, int *et,
double *init, double *trdens, double *dens,
double *alpha, double *beta, double *sca, double *xi) {
@@ -51,40 +51,71 @@
int nttot=0;
- for(int i=0; i<3; i++) {
+ for(int i=0; i<nc[0]; i++) {
nttot += ntimes[i];
- Rprintf("ntimes: %d \n", ntimes[i]);
Rprintf("nttot: %d \n", nttot);
}
// loop over cases
- int tt=0;
+// int tt=0; // main counter should run from 0 to nt-1
for(int cas=0; cas<nc[0]; cas++) {
// compute alpha1 for this case
double sca1=0.0;
+ matrix alpha1(ns[0],1);
+ matrix alphat(ns[0],1);
+ matrix denst(ns[0]);
for(int i=0; i<ns[0]; i++) {
- alpha[tt*ns[0]+i] = init[cas*ns[0]+i]*dens[tt*ns[0]+i];
- sca1 += alpha[tt*ns[0]+i];
+ alpha[(bt[cas]-1)*ns[0]+i] = init[cas*ns[0]+i]*dens[(bt[cas]-1)*ns[0]+i];
+ denst(i+1) = dens[(bt[cas]-1)*ns[0]+i];
+ // compute scale for alpha1
+ sca1 += alpha[(bt[cas]-1)*ns[0]+i];
+// Rprintf("init=%f\n",init[cas*ns[0]+i]);
+// Rprintf("index=%d\n",cas*ns[0]+i);
+// Rprintf("dens=%f\n",dens[(bt[cas]-1)*ns[0]+i]);
+// Rprintf("index=%d\n",(bt[cas]-1)*ns[0]+i);
}
- sca[tt] = sca1;
+ sca[(bt[cas]-1)] = 1/sca1;
+ // scale alpha1
for(int i=0; i<ns[0]; i++) {
- alpha[tt*ns[0]+i] /= sca1;
+ alpha[(bt[cas]-1)*ns[0]+i] /= sca1;
+ alpha1(i+1) = alpha[(bt[cas]-1)*ns[0]+i];
}
-
- // compute scale for alpha1
-
- // scale alpha1
-
+ Rprintf("%d\n",bt[cas]-1);
+ denst.print();
+ alpha1.print();
+ matrix trans(ns[0],ns[0]);
//loop over ntimes[cas]>1
+ if(ntimes[cas]>0) {
+ for(int t=bt[cas]; t<et[cas]; t++) {
+ // get trans and dens values
+ for(int i=0; i<ns[0]; i++) {
+ for(int j=0; j<ns[0]; j++) {
+ trans(i+1,j+1) = trdens[(i*ns[0]+j)*nt[0]+t];
+ }
+ denst(i+1) = dens[t*ns[0]+i];
+ }
+
+ Rprintf("%d\n",t);
+ trans.print();
+// denst.print();
+
+ // compute alphat
+ alphat = had(transpose(trans)*alpha1,denst);
+
+ // compute scale for t
+ sca[t] = 1/(alphat.msum());
+
+ for(int i=0; i<ns[0]; i++) {
+ alphat(i+1) *= sca[t];
+ alpha[t*ns[0]+i] = alphat(i+1);
+ }
+ }
+ }
+
+// Rprintf("%d\n",tt);
- // compute cumulative t for indexing alpha
- if(cas>0) tt += ntimes[cas];
-
- Rprintf("%d\n",tt);
-
} // end cases
-
// (auxiliaries for) gradients from mixing proportions
// matrix *xiC; xiC = new matrix[2];
// mp1[0](1)=3.7;
@@ -93,53 +124,5 @@
}
-
-// //this is the actual loglikelihood function for models with covariates
-// void forwardbackward(double logl, double *init, int *linit) {
-//
-// double crit=0.0;
-//
-// crit=init[0];
-// // This will stop optimization when the loglikelihood is better then crit, useful for
-// // goodness of fit bootstrapping
-// Rprintf("stop crit=%f\n",crit);
-//
-// Rprintf("Starting to compute likelihood.\n");
-//
-// Rprintf("linit: %d \n", linit[0]);
-//
-// logl=0;
-//
-// Rprintf("logl: %f \n", logl);
-//
-// double tmp;
-//
-// for(int i=0; i<3; i++) {
-// logl += init[i];
-// tmp = init[i];
-// Rprintf("init: %f \n", tmp);
-// Rprintf("logl: %f \n", logl);
-// }
-//
-// int nrcomp=models.getNrComp();
-// int ngroups=models.getNrGroups();
-// int *npars = new int[nrcomp];
-// int *nstates = new int[nrcomp];
-// int totalstates=0;
-//
-// for(int cp=0; cp<nrcomp; cp++) {
-// npars[cp]=models.mods[cp].getNPars();
-// nstates[cp]=models.mods[cp].getStates();
-// totalstates += nstates[cp];
-// }
-//
-// // (auxiliaries for) gradients from mixing proportions
-// matrix *mp1; mp1 = new matrix[2];
-// matrix *mpt; mpt = new matrix[ngroups];
-// matrix *mptfinal; mptfinal = new matrix[ngroups];
-
-// }
-
-
} // end extern "C"
\ No newline at end of file
Modified: pkg/depmixS4/src/fb.h
===================================================================
--- pkg/depmixS4/src/fb.h 2011-06-20 08:01:49 UTC (rev 467)
+++ pkg/depmixS4/src/fb.h 2011-06-20 12:42:27 UTC (rev 468)
@@ -15,7 +15,9 @@
#include <Rmath.h>
// criterion for stopping optimization, used in bootstrapping
-void forwardbackward(int *ns, int *nc, int *nt, int *ntimes, double *init, double *trdens, double *dens, double *alpha, double *beta, double *sca, double *xi);
+void forwardbackward(int *ns, int *nc, int *nt, int *ntimes, int *bt, int *et,
+ double *init, double *trdens, double *dens,
+ double *alpha, double *beta, double *sca, double *xi);
// the main function that computes the forward backward vars, and gamma and xi
/* void forwardbackward(double logl, double *init, int *linit); */
More information about the depmix-commits
mailing list