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))) # mxMatrix(name='MP', ncol=2, labels=c('m1Prob', 'm2Prob'), free=F) ) mixRun2 <- mxRun(mixmod2) mixRun2@matrices\$M@values mixRun2@matrices\$D@values mixRun2@objective@result probs2 <- cbind( dpois(eqdat\$V2, mixRun2@matrices\$M@values[,1]), dpois(eqdat\$V2, mixRun2@matrices\$M@values[,2]) ) apply(probs2, 1, which.max) # 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) mixRun3@matrices\$M mixRun3@matrices\$D mixRun3@objective@result