[Depmix-commits] r115 - in trunk: . srcll
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 27 23:26:51 CET 2008
Author: ingmarvisser
Date: 2008-03-27 23:26:51 +0100 (Thu, 27 Mar 2008)
New Revision: 115
Added:
trunk/srcll/
trunk/srcll/.Rhistory
trunk/srcll/filter.R
trunk/srcll/lystig.c
trunk/srcll/lystig.so
trunk/srcll/lystig2.R
trunk/srcll/lystigRfromC.R
Log:
Added srcll directory with tests for C-code for lystig and fb
Added: trunk/srcll/.Rhistory
===================================================================
--- trunk/srcll/.Rhistory (rev 0)
+++ trunk/srcll/.Rhistory 2008-03-27 22:26:51 UTC (rev 115)
@@ -0,0 +1,569 @@
+n
+x
+length(x)
+cfunction
+compileCode
+inline:::compileCode
+?match
+match
+library(dlm)
+?dlm
+## seasonal component for quarterly data#
+dlmModSeas(4, dV = 3.2)
+x <- matrix(runif(6,4,10), nc = 2); x
+dlmModReg(x)
+dlmModReg(x, addInt = FALSE)
+data(NelPlo)#
+### multivariate local level -- seemingly unrelated time series#
+buildSu <- function(x) {#
+ Vsd <- exp(x[1:2])#
+ Vcorr <- tanh(x[3])#
+ V <- Vsd %o% Vsd#
+ V[1,2] <- V[2,1] <- V[1,2] * Vcorr#
+ Wsd <- exp(x[4:5])#
+ Wcorr <- tanh(x[6])#
+ W <- Wsd %o% Wsd#
+ W[1,2] <- W[2,1] <- W[1,2] * Wcorr#
+ return(list(#
+ m0 = rep(0,2),#
+ C0 = 1e7 * diag(2),#
+ FF = diag(2),#
+ GG = diag(2),#
+ V = V,#
+ W = W))#
+}
+head(NelPlo)
+head(NelPlo,30)
+dim(NelPlo)
+NelPlo
+dlmMLE
+dlmLL
+library(rdonlp2)
+library(rdonlp)
+library(Rdonlp)
+library(rdonlp2)
+library(Rdonlp2)
+?Rdonlp2
+?Rdonlp
+?rdonlp
+help(pack="Rdonlp2")
+?donlp2
+library(depmix)
+?depmix
+data(speed)#
+mod <- dmm(nsta=2,itemt=c(1,2)) # gaussian and binary items#
+fit1 <- fitdmm(dat=speed,dmm=mod)
+?donlp2
+library(Rdonlp2)
+y=library(Rdonlp2)
+y
+datasets
+y$datasets
+?datasets
+library(depmix)
+?depmix
+data(speed)#
+mod <- dmm(nsta=2,itemt=c(1,2)) # gaussian and binary items#
+fit1 <- fitdmm(dat
+=speed,dmm=mod)
+summary(fit1)
+library(depmix)
+?depmix
+library(Rhmm)
+?Rhmm
+?HMMSet
+library(RHmm)
+?RHmm
+?HMMSet
+ data(geyser)#
+ obs <- geyser$duration
+head(geyse)
+head(geyser)
+plot(as.ts(geyser$waiting))
+plot(as.ts(geyser$duration))
+plot(as.ts(geyser[1:50,]))
+plot(as.ts(geyser[1:100,]))
+cor(geyser)
+ ResGeyser1 <- HMMFit(obs)
+ ResGeyser2 <- HMMFit(obs, nStates=3, paramBW=list(verbose=1, init="KMEANS"))#
+ # fit a 2 states of a mixture of 3 normal distributions#
+ # for data_mixture#
+ data(data_mixture)#
+ ResMixture <- HMMFit(data_mixture, nStates=2, nMixt=3, dis="MIXTURE")#
+ summary(ResMixture)#
+ # geyser data - 3 states HMM with bivariate normal distribution#
+ ResGeyser<-HMMFit(obs=as.matrix(geyser), nStates=3)#
+ # multiple samples discrete observations#
+ data(weather)#
+ ResDiscrete <- HMMFit(obs=weather, nStates=3, dis="DISCRETE")
+ResDiscrete
+library(nnet)
+?multinom
+multinom
+library(depmixS4)
+?depmix
+ data(speed)#
+ depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()))
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()))
+ mod
+ # create a 2 state model with one continuous and one binary response
+ data(speed)
+ mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()))
+ # print the model, formulae and parameter values
+ mod
+?library
+library(depmixS4)
+?depmix
+?stats4
+library(stats4)
+?stats4
+?statmodel
+library(modeltools)
+?fit
+methods(fit)
+getAnywhere(fit)
+fit[depmixS4]
+fit[1]
+[2]
+modeltools:::fit
+depmixS4:::fit
+depmix
+?depmix
+?depmix.fit
+fit
+ data(speed)
+ # 2-state model on the RTs of the speed data with random
+ # starting values for the transition pars (without those EM does not get off the ground)
+ set.seed(1)
+ mod <- depmix(rt~1,data=speed,nstates=2,trstart=runif(4))
+ # fit the model
+ mod1 <- fit(mod)
+ mod1 # to see the logLik and optimization information
+ # to see the parameters
+ summary(mod1)
+mod
+library(modeltools)
+library(depmixS4)
+?depmix.fit
+ls()
+?depmix-class
+?"depmix-class"
+library(depmixS4)
+?"depmix-class"
+?depmix.fit
+library(depmixS4)
+?depmix.fit
+data(speed)#
+# 2-state model on the RTs of the speed data with random #
+# starting values for the transition pars (without those EM does not get off the ground)#
+set.seed(1)#
+mod <- depmix(rt~1,data=speed,nstates=2,trstart=runif(4))#
+# fit the model#
+mod1 <- fit(mod)
+mod1
+summary(mod1)
+?depmix
+?llratio
+library(depmixS4)
+?llratio
+?depmix
+library(depmixS4)
+?llratio
+library(depmixS4)
+?llratio
+library(depmixS4)
+?llratio
+?balance
+posterior
+?balance
+library(flexmix)
+?flexmix
+?balance
+?flexmix
+?response-class
+?"response-class"
+?speed
+data(speed)
+head(speed)
+head(speed,24)
+max(speed$Pacc)
+?"response-classes"
+?"response-class"
+?glm
+library(depmixS4)
+?balance
+?"depmix-class"
+?em
+?"depmix.fi"
+?"depmix.fit"
+?depmix
+mod
+
+# to see the ordering of parameters to use in setpars
+mod <- setpars(mod, value=1:npar(mod))
+mod
+# create a 2 state model with one continuous and one binary response
+data(speed)
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()))
+# print the model, formulae and parameter values
+mod
+
+# to see the ordering of parameters to use in setpars
+mod <- setpars(mod, value=1:npar(mod))
+mod
+mod <- setpars(mod, getpars(mod,which="fixed"))
+mod
+?llratio
+?"response-class"
+?GLMresponse-class
+?"GLMresponse-class"
+?glm
+?response
+library(depmixS4)
+
+x <- rnorm(100)
+xd <- data.frame(x,1)
+
+mod <- depmix(x~1,ns=2,nt=100,trst=runif(4))
+
+viterbi(mod)
+fm <- fit(mod)
+
+viterbi(fm)
+data(speed)
+rt <- speed$rt
+
+mod <- depmix(rt~1,ns=2,nt=439,trst=runif(4))
+
+fm <- fit(mod)
+
+viterbi(fm)
+viterbi
+library(depmixS4)
+data(speed)
+rt <- speed$rt
+
+mod <- depmix(rt~1,ns=2,nt=439,trst=runif(4))
+
+fm <- fit(mod)
+
+viterbi(fm)
+viterbi
+library(depmixS4)
+data(speed)
+rt <- speed$rt
+
+mod <- depmix(rt~1,ns=2,nt=439,trst=runif(4))
+
+fm <- fit(mod)
+
+viterbi(fm)
+plot(as.ts(viterbi(fm))
+)
+plot(as.ts(viterbi(fm)))
+trstart=c(0.899,0.101,0.084,0.916)
+instart=c(0.5,0.5)
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()),trstart=trstart)
+# respstart=resp,trstart=trstart,instart=instart)
+
+logLik(mod)
+
+mod1 <- fit(mod)
+ll <- logLik(mod1)
+
+
+#
+# Test optimization using Rdonlp2
+#
+
+trstart=c(0.899,0.101,0,0.01,0.084,0.916,0,0)
+instart=c(0.5,0.5)
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,family=list(gaussian(),multinomial()),
+ respstart=resp,trstart=trstart,instart=instart)
+
+logLik(mod)
+
+mod1 <- fit(mod)
+ll <- logLik(mod1)
+post <- cbind(viterbi(mod1),speed$Pacc)
+plot(as.ts(post))
+cor(post)
+mod at ntimes
+mod1 at ntimes
+plot(rnorm(100))#
+par(fig=c(0, 1/2, 0, 1/2), new=T)#
+plot(seq(-2,2,length=300),dnorm(seq(-2,2,length=300)),type="l", axes =#
+F, xlab="", ylab="")
+A <- matrix(1:4,2,2)
+B <- matrix(1:10,10)
+init <- matrix(1:2,1)
+A
+B
+init
+ct <- A*B
+apply(B,1,prod,A)
+apply(B,1,,A)
+A*B[1,]
+apply(B,1,*,A)
+apply(B,1,"*",A)
+?array
+ct <- array(apply(B,1,"*",A),c(2,2,10))
+ct
+library(depmixS4)
+
+data(speed)
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()))
+# print the model, formulae and parameter values
+mod
+
+x=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+x
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/srcll")
+source("lystig2.R")
+x1=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+
+source("lystig2.R")
+
+x2=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+
+all.equal(x1,x2)
+x1=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+
+source("lystig2.R")
+
+x2=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+
+all.equal(x1,x2)
+x1=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+x1b=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),439)
+
+
+source("lystig2.R")
+
+x2=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+x2b=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),439)
+
+all.equal(x1,x2)
+all.equal(x1b,x2b)
+A <- matrix(1:4,2,2)
+B <- matrix(1:10,5)
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,10))
+ct
+B
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,10))
+A
+B
+ct
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,5))
+ct
+A
+B
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,5),byrow=T)
+?array
+methods(summary)
+showMethods(summary)
+A[1,]*B[1,]
+A[2,]*B[1,]
+apply(B,1,"*",A),c(2,2,5)
+apply(B,1,"*",A)
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,5,2)
+)
+ct
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,5))
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,5))
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct[1,]
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(2,2,5))
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct[1,,]
+
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct[1,,]
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+A[1,]*B[1,]
+A[2,]*B[1,]
+A <- matrix(1:4,2,2)
+
+B <- matrix(1:10,5)
+
+init <- matrix(1:2,1)
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+ct <- array(apply(B,1,"*",A),c(2,5,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+ct <- array(apply(B,1,"*",A),c(2,2,5))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+ct <- array(apply(A,1,"*",B),c(2,2,5))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+ct <- array(apply(A,1,"*",B),c(2,5,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+ct <- array(apply(A,1,"*",B),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+A[1,]*B[1,]
+A[2,]*B[1,]
+ct <- array(apply(A,1,"*",B),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+
+
+A[1,]*B[1,]
+A[2,]*B[1,]
+t(ct)
+ct <- array(apply(A,1,"*",B),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+
+
+A[1,]*B[1,]
+A[2,]*B[1,]
+
+A[1,]*B[2,]
+A[2,]*B[2,]
+ct[2,,]
+ct <- array(apply(t(A),1,"*",B),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+
+
+A[1,]*B[1,]
+A[2,]*B[1,]
+ct <- array(apply(A,1,"*",B),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+
+
+A[1,]*B[1,]
+A[2,]*B[1,]
+sessionInfo()
+A <- array(1:20,c(5,2,2))
+
+B <- matrix(1:10,5)
+A
+B
+A <- array(1:20,c(2,2,5))
+A
+B <- matrix(1:10,5)
+B
+rbind(A[1,]*B[1,],
+A[2,]*B[1,])
+
+rbind(A[1,,1]*B[1,],
+A[2,,1]*B[1,])
+rbind(A[1,,2]*B[2,],
+A[2,,2]*B[2,])
+rbind(A[1,,3]*B[2,],
+A[2,,3]*B[2,])
+apply(B,1,"*",A)
+c(sort(apply(B,1,"*",A)))
+ct <- numeric(0)
+for(i in 1:5){
+ct <- c(ct,(A[1,,i]*B[1,],
+A[2,,i]*B[1,])
+}
+ct <- numeric(0)
+ct
+ct <- numeric(0)
+for(i in 1:5){
+ct <- c(ct,A[1,,i]*B[1,],A[2,,i]*B[1,])
+}
+ct
Added: trunk/srcll/filter.R
===================================================================
--- trunk/srcll/filter.R (rev 0)
+++ trunk/srcll/filter.R 2008-03-27 22:26:51 UTC (rev 115)
@@ -0,0 +1,36 @@
+#
+# Ingmar Visser
+#
+# lystig filter set up
+#
+
+init <- matrix(1:2,1)
+
+
+A <- array(1:20,c(2,2,5))
+
+B <- matrix(1:10,5)
+
+
+rbind(A[1,,1]*B[1,],
+A[2,,1]*B[1,])
+
+rbind(A[1,,2]*B[2,],
+A[2,,2]*B[2,])
+
+rbind(A[1,,3]*B[2,],
+A[2,,3]*B[2,])
+
+
+
+ct <- array(apply(B,1,"*",A),c(5,2,2))
+
+ct[1,,]
+ct[,1,]
+ct[,,1]
+
+ct <- numeric(0)
+for(i in 1:5){
+ct <- c(ct,A[1,,i]*B[1,],A[2,,i]*B[1,])
+}
+
Added: trunk/srcll/lystig.c
===================================================================
--- trunk/srcll/lystig.c (rev 0)
+++ trunk/srcll/lystig.c 2008-03-27 22:26:51 UTC (rev 115)
@@ -0,0 +1,79 @@
+/*
+moving lystig into C code
+*/
+
+#include <R.h>
+#include <Rinternals.h>
+
+/*
+ * Log likelihood computation according to Lystig & Hughes (2002). This
+ * is very similar to the Forward part of the Forward-Backward algorithm
+ * but admits of easy computation of the gradients of parameters and the
+ * observed information. This version of the routine only computes the
+ * likelihood though.
+ *
+ * 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 = K*K matrix with transition probabilities, from column to row !!!!!!!
+ * change to T*K*K
+ *
+ * B = T*K matrix with elements ab_{ij} = P(y_i|s_j)
+ * init = K vector with initial probabilities !!!OR!!! K*length(ntimes) matrix with initial probs per case
+ */
+
+/*
+ * RETURNS: 'sca'le factors, recurrent variables 'phi', loglikelihood
+ */
+
+
+
+
+SEXP mkans(double x) {
+ SEXP ans;
+ PROTECT(ans = allocVector(REALSXP, 1));
+ REAL(ans)[0] = x;
+ UNPROTECT(1);
+ return ans;
+}
+
+/* , SEXP A, SEXP B, SEXP, SEXP ntimes, SEXP stationary */
+
+SEXP lystig(SEXP init, SEXP ntimes) {
+
+ int nt=INTEGER(ntimes)[0];
+
+ Rprintf("nt: %d \n", nt);
+
+ double loglike=0.0;
+
+ Rprintf("Init ans: %f \n",loglike);
+
+ for(int i=0; i<nt; i++) {
+ loglike+=1;
+ }
+
+ Rprintf("End ans: %f \n",loglike);
+
+ return(mkans(loglike));
+
+}
+
+
+
+/*
+ * for(case in 1:lt) { # multiple cases
+ * phi[bt[case],] <- init[case,]*B[bt[case],] # initialize
+ * sca[bt[case]] <- 1/sum(phi[bt[case],])
+ * if(ntimes[case]>1) {
+ * for(i in (bt[case]+1):et[case]) {
+ * if(stationary) phi[i,] <- (A[1,,]%*%phi[i-1,])*B[i,]
+ * else phi[i,] <- (A[i-1,,]%*%phi[i-1,])*B[i,]
+ * phi[i,] <- sca[i-1]*phi[i,]
+ * sca[i] <- 1/sum(phi[i,])
+ * }
+ * }
+ * }
+ *
+ * logLike=-sum(log(sca))
+ * return(list(phi=phi,sca=sca,logLike=logLike))
+ */
Property changes on: trunk/srcll/lystig.c
___________________________________________________________________
Name: svn:executable
+ *
Added: trunk/srcll/lystig.so
===================================================================
(Binary files differ)
Property changes on: trunk/srcll/lystig.so
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mime-type
+ application/octet-stream
Added: trunk/srcll/lystig2.R
===================================================================
--- trunk/srcll/lystig2.R (rev 0)
+++ trunk/srcll/lystig2.R 2008-03-27 22:26:51 UTC (rev 115)
@@ -0,0 +1,61 @@
+#
+# Ingmar Visser
+#
+# LYSTIG algoritme voor de loglikelihood, 23-3-2008
+#
+
+lystig <- function(init,A,B,ntimes=NULL,stationary=TRUE) {
+
+ # Log likelihood computation according to Lystig & Hughes (2002). This
+ # is very similar to the Forward part of the Forward-Backward algorithm
+ # but admits of easy computation of the gradients of parameters and the
+ # observed information. This version of the routine only computes the
+ # likelihood though.
+
+ # NOTE THE CHANGE IN FROM ROW TO COLUMN SUCH THAT TRANSPOSING A IS NOT
+ # NECCESSARY ANYMORE
+
+ # K is the number of states of the model
+
+ # A = K*K matrix if stationary=TRUE with transition probabilities, from column to row!!!!!
+ # A = T*K*K arrary if stationary=FALSE
+
+ # B = T*K matrix with elements b_{tj} = P(y_t|s_j)
+ # init = K*length(ntimes) matrix with initial probs per case
+ # Returns: 'sca'le factors, recurrent variables 'phi', loglikelihood
+
+ nt <- nrow(B)
+ ns <- ncol(init)
+
+ if(!is.null(ntimes)) {
+ ntimes <- nt
+ }
+
+ phi <- matrix(ncol=ns,nrow=nt)
+ sca <- vector(length=nt)
+
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ ll <- 0
+
+ for(case in 1:lt) { # multiple cases
+ phi[bt[case],] <- init[case,]*B[bt[case],] # initialize
+ sca[bt[case]] <- 1/sum(phi[bt[case],])
+ if(ntimes[case]>1) {
+ for(i in (bt[case]+1):et[case]) {
+ if(stationary) phi[i,] <- (A[1,,]%*%phi[i-1,])*B[i,]
+ else phi[i,] <- (A[i-1,,]%*%phi[i-1,])*B[i,]
+ phi[i,] <- sca[i-1]*phi[i,]
+ sca[i] <- 1/sum(phi[i,])
+ }
+ }
+ }
+
+ logLike=-sum(log(sca))
+
+ return(list(phi=phi,sca=sca,logLike=logLike))
+
+}
+
Added: trunk/srcll/lystigRfromC.R
===================================================================
--- trunk/srcll/lystigRfromC.R (rev 0)
+++ trunk/srcll/lystigRfromC.R 2008-03-27 22:26:51 UTC (rev 115)
@@ -0,0 +1,38 @@
+
+cd /Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/srcll
+
+R CMD SHLIB lystig.c
+
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/srcll")
+
+dyn.load("lystig.so")
+
+#
+# also look at filter!!!!!
+#
+
+library(depmixS4)
+
+data(speed)
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()))
+# print the model, formulae and parameter values
+mod
+
+x1=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+x1b=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),439)
+
+
+source("lystig2.R")
+
+x2=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),mod at ntimes)
+x2b=lystig(mod at init,mod at trDens,apply(mod at dens,c(1,3),prod),439)
+
+all.equal(x1,x2)
+all.equal(x1b,x2b)
+
+
+
+.Call("lystig",as.integer(mod at ntimes))
+
+dyn.unload("lystig.so")
+
More information about the depmix-commits
mailing list