[Depmix-commits] r512 - pkg/depmixS4/R

[noreply at r-forge.r-project.org]

Author: maarten
Date: 2012-03-29 12:02:43 +0200 (Thu, 29 Mar 2012)
New Revision: 512

Modified:
   pkg/depmixS4/R/EM.R
Log:
- random gamma starting values in EM now uses a dirichlet distribution
- TODO: specify dispersion as a parameter in em.options()


Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R	2012-03-16 13:17:11 UTC (rev 511)
+++ pkg/depmixS4/R/EM.R	2012-03-29 10:02:43 UTC (rev 512)
@@ -2,6 +2,15 @@
 # Maarten Speekenbrink 23-3-2008
 # 
 
+rdirichlet <- function(n, alpha) {
+  # taken from gtools...
+    l <- length(alpha)
+    x <- matrix(rgamma(l * n, alpha), ncol = l, byrow = TRUE)
+    sm <- x %*% rep(1, l)
+    x/as.vector(sm)
+}
+
+
 em <- function(object,...) {
 	if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
 	call <- match.call()
@@ -31,8 +40,10 @@
 	if(random.start) {
 				
 		nr <- sum(ntimes(object))
-		gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nrow=nr,ncol=ns)
-		gamma <- gamma/rowSums(gamma)
+		#gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nrow=nr,ncol=ns)
+		#gamma <- matrix(rbeta(nr*ns,alpha=.1,beta=.1),nrow=nr,ncol=ns)
+		gamma <- rdirichlet(nr,alpha=rep(.1,ns))
+		#gamma <- gamma/rowSums(gamma)
 		LL <- -1e10
 		
 		for(i in 1:ns) {
@@ -143,8 +154,9 @@
 	if(random.start) {
 				
 		nr <- sum(ntimes(object))
-		gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nrow=nr,ncol=ns)
-		gamma <- gamma/rowSums(gamma)
+		#gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nrow=nr,ncol=ns)
+		#gamma <- gamma/rowSums(gamma)
+		gamma <- rdirichlet(nr,alpha=rep(.1,ns))
 		LL <- -1e10
 		
 		for(i in 1:ns) {