[Depmix-commits] r465 - in pkg/depmixS4: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jun 19 23:07:08 CEST 2011
Author: ingmarvisser
Date: 2011-06-19 23:07:06 +0200 (Sun, 19 Jun 2011)
New Revision: 465
Modified:
pkg/depmixS4/R/allGenerics.R
pkg/depmixS4/R/fb.R
pkg/depmixS4/R/forwardbackward.R
pkg/depmixS4/src/fb.cc
pkg/depmixS4/src/fb.h
Log:
Initial set up for C code of forward backward algorithm.
Modified: pkg/depmixS4/R/allGenerics.R
===================================================================
--- pkg/depmixS4/R/allGenerics.R 2011-06-17 13:25:47 UTC (rev 464)
+++ pkg/depmixS4/R/allGenerics.R 2011-06-19 21:07:06 UTC (rev 465)
@@ -6,7 +6,7 @@
# 17-6-2011: added dynamic lib statement to include the C code
# version of forward backward routine
-.First.lib <- function(lib, pkg) {
+.onLoad <- function(lib, pkg) {
require(stats)
require(methods)
require(MASS)
@@ -16,7 +16,7 @@
library.dynam("depmixS4", pkg, lib)
}
-.Last.lib <- function(libpath) {
+.onUnLoad <- function(libpath) {
library.dynam.unload("depmixS4",libpath)
}
Modified: pkg/depmixS4/R/fb.R
===================================================================
--- pkg/depmixS4/R/fb.R 2011-06-17 13:25:47 UTC (rev 464)
+++ pkg/depmixS4/R/fb.R 2011-06-19 21:07:06 UTC (rev 465)
@@ -4,7 +4,7 @@
# FORWARD-BACKWARD algoritme, 23-3-2008
#
-fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE) {
+fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE,useC=TRUE) {
# Forward-Backward algorithm (used in Baum-Welch)
# Returns alpha, beta, and full data likelihood
@@ -19,61 +19,82 @@
# NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i) !!!NOTE the order of i and j!!!
+
+
+
B <- apply(B,c(1,3),prod)
nt <- nrow(B)
ns <- ncol(init)
- alpha <- matrix(ncol=ns,nrow=nt)
- beta <- matrix(ncol=ns,nrow=nt)
- sca <- vector(length=nt)
- xi <- array(dim=c(nt,ns,ns))
-
+ 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)
- et <- cumsum(ntimes)
- bt <- c(1,et[-lt]+1)
- for(case in 1:lt) {
- alpha[bt[case],] <- init[case,]*B[bt[case],] # initialize
- sca[bt[case]] <- 1/sum(alpha[bt[case],])
- alpha[bt[case],] <- alpha[bt[case],]*sca[bt[case]]
+ if(useC) {
+ 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")
- if(ntimes[case]>1) {
- for(i in bt[case]:(et[case]-1)) {
- if(stationary) alpha[i+1,] <- (A[1,,]%*%alpha[i,])*B[i+1,]
- else alpha[i+1,] <- (A[i,,]%*%alpha[i,])*B[i+1,]
- sca[i+1] <- 1/sum(alpha[i+1,])
- alpha[i+1,] <- sca[i+1]*alpha[i+1,]
- }
- }
+ } else {
- beta[et[case],] <- 1*sca[et[case]] # initialize
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
- if(ntimes[case]>1) {
- for(i in (et[case]-1):bt[case]) {
- if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
- else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
+ for(case in 1:lt) {
+ alpha[bt[case],] <- init[case,]*B[bt[case],] # initialize
+ sca[bt[case]] <- 1/sum(alpha[bt[case],])
+ alpha[bt[case],] <- alpha[bt[case],]*sca[bt[case]]
+
+ if(ntimes[case]>1) {
+ for(i in bt[case]:(et[case]-1)) {
+ if(stationary) alpha[i+1,] <- (A[1,,]%*%alpha[i,])*B[i+1,]
+ else alpha[i+1,] <- (A[i,,]%*%alpha[i,])*B[i+1,]
+ sca[i+1] <- 1/sum(alpha[i+1,])
+ alpha[i+1,] <- sca[i+1]*alpha[i+1,]
+ }
}
- for(i in bt[case]:(et[case]-1)) {
- if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
- else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
+ beta[et[case],] <- 1*sca[et[case]] # initialize
+
+ if(ntimes[case]>1) {
+ for(i in (et[case]-1):bt[case]) {
+ if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
+ else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
+ }
+
+ for(i in bt[case]:(et[case]-1)) {
+ if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
+ else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
+ }
}
+
}
+
+ gamma <- alpha*beta/sca
+ like <- -sum(log(sca))
+ if(return.all) {
+ res <- list(alpha=alpha,beta=beta,gamma=gamma,xi=xi,sca=sca,logLike=like)
+ } else {
+ res <- list(gamma=gamma,xi=xi,logLike=like)
+ }
}
-
- gamma <- alpha*beta/sca
- like <- -sum(log(sca))
-
- if(return.all) {
- res <- list(alpha=alpha,beta=beta,gamma=gamma,xi=xi,sca=sca,logLike=like)
- } else {
- res <- list(gamma=gamma,xi=xi,logLike=like)
- }
-
+
res
}
Modified: pkg/depmixS4/R/forwardbackward.R
===================================================================
--- pkg/depmixS4/R/forwardbackward.R 2011-06-17 13:25:47 UTC (rev 464)
+++ pkg/depmixS4/R/forwardbackward.R 2011-06-19 21:07:06 UTC (rev 465)
@@ -5,9 +5,9 @@
#
setMethod("forwardbackward","depmix",
- function(object, return.all=TRUE, ...) {
+ function(object, return.all=TRUE, useC=TRUE, ...) {
fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),
- stationary=object at stationary,return.all=return.all)
+ stationary=object at stationary,return.all=return.all,useC=useC)
}
)
Modified: pkg/depmixS4/src/fb.cc
===================================================================
--- pkg/depmixS4/src/fb.cc 2011-06-17 13:25:47 UTC (rev 464)
+++ pkg/depmixS4/src/fb.cc 2011-06-19 21:07:06 UTC (rev 465)
@@ -1,26 +1,124 @@
#include "fb.h"
+#include "matrix.h"
+
extern "C" {
+
+// # Forward-Backward algorithm (used in Baum-Welch)
+// # Returns alpha, beta, and full data likelihood
+//
+// # NOTE THE CHANGE IN FROM ROW TO COLUMN SUCH THAT TRANSPOSING A IS NOT NECCESSARY ANYMORE
+// # IN COMPUTING ALPHA AND BETA BUT IS NOW NECCESSARY IN COMPUTING XI
+// # A = T*K*K matrix with transition probabilities, from row to column!!!!!!!
+// # B = T*K matrix with elements ab_{ij} = P(y_i|s_j)
+// # init = K vector with initial probabilities
+//
+// # NOTE: to prevent underflow, alpha and beta are scaled, using sca
+//
+// # NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i) !!!NOTE the order of i and j!!!
-//this is the actual loglikelihood function for models with covariates
-double forwardbackward(double *init, int *linit) {
+// fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE) {
+//
+// fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),
+// stationary=object at stationary,return.all=return.all)
- Rprintf("Starting to compute likelihood.\n");
-
- int ltin;
-
- ltin = linit[0];
-
-// Rprintf(ltin);
-
- double res;
- res=0;
-
- for(int i=0; i<ltin; i++) {
- res += init[i];
+
+// inputs are:
+// a) ns: the number of states
+// b) nc: the number of cases
+// c) nt: the number of rows of data
+// d) ntimes: rows of data of individual cases
+// 1) init: ns vector or ns by nc matrix
+// 2) trdens: ns by ns matrix or nt by ns by ns array
+// 3) dens: nt by ns matrix
+//
+// outputs are:
+// 1) alpha: nt by ns matrix
+// 2) beta: nt by ns matrix
+// 3) xi: nt by ns by ns array
+// 4) sca: nt vector
+
+// gamma is computed as alpha*beta/sca in R (no loop needed)
+
+
+void forwardbackward(int *ns, int *nc, int *nt, int *ntimes, double *init, double *trdens, double *dens, double *alpha, double *beta, double *sca, double *xi) {
+
+ Rprintf("ns=%d\n",ns[0]);
+ Rprintf("nc=%d\n",nc[0]);
+ Rprintf("nt=%d\n",nt[0]);
+
+ int nttot=0;
+
+ for(int i=0; i<3; i++) {
+ nttot += ntimes[i];
+ Rprintf("ntimes: %d \n", ntimes[i]);
+ Rprintf("nttot: %d \n", nttot);
+ }
+
+ // loop over cases
+ int tt=0;
+ for(int cas=0; cas<nc[0]; cas++) {
+ // compute alpha1 for this case
+ double sca1=0.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];
+ }
+ sca[tt] = sca1;
+ for(int i=0; i<ns[0]; i++) {
+ alpha[tt*ns[0]+i] /= sca1;
+ }
+
+ // compute scale for alpha1
+
+ // scale alpha1
+
+ //loop over ntimes[cas]>1
+
+ // 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;
+// mp1[1](1)=1.6;
+// mp1[0].print();
+
}
+
-
+// //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];
@@ -37,10 +135,9 @@
// matrix *mp1; mp1 = new matrix[2];
// matrix *mpt; mpt = new matrix[ngroups];
// matrix *mptfinal; mptfinal = new matrix[ngroups];
-
- return(res);
-
-}
+
+// }
+
} // end extern "C"
\ No newline at end of file
Modified: pkg/depmixS4/src/fb.h
===================================================================
--- pkg/depmixS4/src/fb.h 2011-06-17 13:25:47 UTC (rev 464)
+++ pkg/depmixS4/src/fb.h 2011-06-19 21:07:06 UTC (rev 465)
@@ -1,17 +1,25 @@
+#ifndef FB
+#define FB 1
+
#include <stdio.h>
#include <stdlib.h>
/* #include <fstream.h> */
-#include "matrix.h"
-
+
+/* #include "matrix.h" */
+
extern "C" {
-
-#include <R.h>
+
+#include <R.h>
#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);
+
// the main function that computes the forward backward vars, and gamma and xi
-double forwardbackward(double *init, int *linit);
-
+/* void forwardbackward(double logl, double *init, int *linit); */
+
} //end extern "C"
+#endif
\ No newline at end of file
More information about the depmix-commits
mailing list