[Picante-commits] r126 - branches/gsoc/R

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

Author: mrhelmus
Date: 2008-06-29 16:59:36 +0200 (Sun, 29 Jun 2008)
New Revision: 126

Modified:
   branches/gsoc/R/pblm.r
Log:
added maxit to function

Modified: branches/gsoc/R/pblm.r
===================================================================
--- branches/gsoc/R/pblm.r	2008-06-29 09:04:34 UTC (rev 125)
+++ branches/gsoc/R/pblm.r	2008-06-29 14:59:36 UTC (rev 126)
@@ -1,4 +1,4 @@
-pblm<-function(assocs,tree1=NULL,tree2=NULL,covars1=NULL,covars2=NULL,bootstrap=FALSE){
+pblm<-function(assocs,tree1=NULL,tree2=NULL,covars1=NULL,covars2=NULL,bootstrap=FALSE,maxit=10000){
 
 
   # Make a vector of associations
@@ -141,8 +141,9 @@
     V1<-as.matrix(V1)
     V2<-as.matrix(V2)
     
-    V1<-V1/det(V1)^(1/nspp1)   # model of coevolution
-  	V2<-V2/det(V2)^(1/nspp2)
+    
+    V1<-V1/det(V1)^(1/nspp1)   # scale covariance matrices (this reduces numerical problems caused by
+  	V2<-V2/det(V2)^(1/nspp2)   # determinants going to infinity or zero)
   	V<-kronecker(V2,V1)  
     invV<-qr.solve(V)
     
@@ -184,7 +185,7 @@
       t(E)%*%invV%*%E/(nassocs-1)
     }
     # estimate d1 and d2 via Nelder-Mead method same as fminsearch in Matlab, by minimizing MSE
-    est<-optim(pstart,pegls,control=list(maxit=10000))        
+    est<-optim(pstart,pegls,control=list(maxit=maxit))        
     MSEFull<-est$value
   	d1<-abs(est$par[1])
   	d2<-abs(est$par[2])