R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> 
> eqdat <- read.table('~/Documents/Projects/HiddenMarkovModelOMx2010/ZucchiniMacDonald2009HMMwithR/dataEarthquakes.txt')
> library(OpenMx)
> 
> #----------------------------------------------------
> # Independent Mixture Distribution Models
> 
> # 2 parameter model
> mixmod2 <- mxModel(
+ 	name='Model to Calculate Poisson Mixture with 2 Distributions',
+ 	mxData(eqdat, 'raw'),
+ 	#
+ 	# Make a row matrix, D, that is the initial distribution: the baseline probability of each latent distribution
+ 	mxMatrix(
+ 		name='D',
+ 		values=.5,
+ 		lbound=0,
+ 		ubound=1,
+ 		free=T,
+ 		nrow=1,
+ 		ncol=2),
+ 	#
+ 	# Constrain D to be a probability vector, i.e. the probability of SOMETHING happening is 1
+ 	# The sum of all the entries should be 1
+ 	mxAlgebra(name='DSUM', sum(D[1,])),
+ 	mxConstraint(DSUM==1),
+ 	#
+ 	# Make a matrix, M, that has the means to be estimated for the latent distributions
+ 	mxMatrix(
+ 		name='M',
+ 		nrow=1,
+ 		ncol=2,
+ 		values=c(16, 24),
+ 		free=T,
+ 		labels=c('m1', 'm2')),
+ 	#
+ 	# Calculate the probabilty of the data given a Poisson distribution with mean m1
+ 	mxAlgebra(name='m1Prob', m1^data.V2 * exp(-m1) / prod(data.V2:1)),
+ 	#
+ 	# Calculate the probabilty of the data given a Poisson distribution with mean m2
+ 	mxAlgebra(name='m2Prob', m2^data.V2 * exp(-m2) / prod(data.V2:1)),
+ 	#
+ 	# Combine in two columns the probabilities calculated above
+ 	mxAlgebra(name='MProb', cbind(m1Prob, m2Prob)),
+ 	#
+ 	# Make the weighted combination of the probs, based on the initial distribution (base rates)
+ 	mxAlgebra(name='WC', MProb%*%t(D)), # Weighted Combination
+ 	mxRowObjective('WC', 'Res', 'Red'),
+ 	mxAlgebra(name='Red', -2*sum(log(Res)))
+ )
> 
> mixRun2 <- mxRun(mixmod2)
Running Model to Calculate Poisson Mixture with 2 Distributions 
> mixRun2@matrices$M
FullMatrix 'M' 

@labels
     [,1] [,2]
[1,] "m1" "m2"

@values
         [,1]     [,2]
[1,] 15.77711 26.83989

@free
     [,1] [,2]
[1,] TRUE TRUE

@lbound: No lower bounds assigned.

@ubound: No upper bounds assigned.

> mixRun2@matrices$D
FullMatrix 'D' 

@labels: No labels assigned.

@values
          [,1]      [,2]
[1,] 0.6757257 0.3242743

@free
     [,1] [,2]
[1,] TRUE TRUE

@lbound
     [,1] [,2]
[1,]    0    0

@ubound
     [,1] [,2]
[1,]    1    1

> 
> mixRun2@objective@result
         [,1]
[1,] 720.7381
> 
> probs2 <- cbind(
+ 	dpois(eqdat$V2, mixRun2@matrices$M@values[,1]),
+ 	dpois(eqdat$V2, mixRun2@matrices$M@values[,2])
+ 	)
> apply(probs2, 1, which.max)
  [1] 1 1 1 1 1 2 2 2 1 2 2 2 2 2 2 1 2 2 2 1 1 1 1 2 1 1 1 1 2 1 1 2 1 1 2 2 2
 [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 1 2 1 1 2 1 1 1 1 2 1 1 2 2 2 2 1 1
 [75] 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> 
> 
> # 3 parameter model
> mixmod3 <- mxModel(
+ 	name='Model to Calculate Poisson Mixture with 3 Distributions',
+ 	mxData(eqdat, 'raw'),
+ 	mxMatrix(
+ 		name='D',
+ 		values=1/3,
+ 		lbound=0,
+ 		ubound=1,
+ 		free=T,
+ 		nrow=1,
+ 		ncol=3),
+ 	mxAlgebra(name='DSUM', sum(D[1,])),
+ 	mxConstraint(DSUM==1),
+ 	mxMatrix(
+ 		name='M',
+ 		nrow=1,
+ 		ncol=3,
+ 		values=c(16, 24, 28),
+ 		free=T,
+ 		labels=c('m1', 'm2', 'm3')),
+ 	mxAlgebra(name='m1Prob', m1^data.V2 * exp(-m1) / prod(data.V2:1)),
+ 	mxAlgebra(name='m2Prob', m2^data.V2 * exp(-m2) / prod(data.V2:1)),
+ 	mxAlgebra(name='m3Prob', m3^data.V2 * exp(-m3) / prod(data.V2:1)),
+ 	mxAlgebra(name='MProb', cbind(m1Prob, m2Prob, m3Prob)),
+ 	mxAlgebra(name='WC', MProb%*%t(D)), # Weighted Combination 
+ 	mxRowObjective('WC', 'Res', 'Red'),
+ 	mxAlgebra(name='Red', -2*sum(log(Res)))
+ )
> 
> mixRun3 <- mxRun(mixmod3)
Running Model to Calculate Poisson Mixture with 3 Distributions 
> 
> mixRun3@matrices$M
FullMatrix 'M' 

@labels
     [,1] [,2] [,3]
[1,] "m1" "m2" "m3"

@values
         [,1]     [,2]     [,3]
[1,] 12.73578 19.78522 31.62954

@free
     [,1] [,2] [,3]
[1,] TRUE TRUE TRUE

@lbound: No lower bounds assigned.

@ubound: No upper bounds assigned.

> mixRun3@matrices$D
FullMatrix 'D' 

@labels: No labels assigned.

@values
          [,1]      [,2]      [,3]
[1,] 0.2775381 0.5928017 0.1296602

@free
     [,1] [,2] [,3]
[1,] TRUE TRUE TRUE

@lbound
     [,1] [,2] [,3]
[1,]    0    0    0

@ubound
     [,1] [,2] [,3]
[1,]    1    1    1

> 
> mixRun3@objective@result
         [,1]
[1,] 713.6979
>