Replace compare_slope_generic_predict_f.R for handling situations with only one predictor

R/compare_slope_generic_predict_f.R

@@ -20,14 +20,18 @@ compare_slope_generic_predict_f = function(list_global, model1 = FALSE)
   		as.numeric(interpolate.pred(rbind(c(M1$fit[length(M1$fit)],M1$pred)),0.01))[-1])
   }
   
-  apply(Ztrunc[-(1:2),],1,function(u) mean(abs(Ztrunc[1,]-u)))
-  apply(Ztrunc[-(1:2),],1,function(u) sd(Ztrunc[1,]-u))
+  #browser()
+  
+  #apply(Ztrunc[-(1:2),],1,function(u) mean(abs(Ztrunc[1,]-u)))
+  #apply(Ztrunc[-(1:2),],1,function(u) sd(Ztrunc[1,]-u))
   
   
   series0=1
   seriescompet=3:nrow(Ztrunc)
   Ncompet=length(seriescompet)
-  Ncompet
+  #Ncompet
+  
+  
   
   ## contrib - poisson
   if(mixture=="dpois"){
@@ -93,7 +97,7 @@ compare_slope_generic_predict_f = function(list_global, model1 = FALSE)
   if(mixture=="dnbinom"){
     loglik=function(Z0,Z,p,sigma,penal.lambda){
       D0=Z0[-1]-Z0[-length(Z0)]
-      D=Z[,-1]-Z[,-length(Z0)]
+      D=rbind(Z[,-1])-rbind(Z[,-length(Z0)])
       pr=sapply(1:nrow(Z),function(i) p[i]*dnbinom(D0,mu=D[i,],size=sigma[i]))
       colNA=(1:ncol(pr))[is.na(colMeans(pr))]
       for(i in colNA){
@@ -128,28 +132,47 @@ compare_slope_generic_predict_f = function(list_global, model1 = FALSE)
       pred2=pred2[,colMeans(!is.na(pred2))>=0.5]
       return(list(pred=pred,predmono=pred2))
     }
-    p.init=1/apply(Ztrunc[-(1:2),],1,function(u) mean((Ztrunc[1,]-u)^2,na.rm=TRUE))
+    
+    
+    p.init=1/apply(rbind(Ztrunc[-(1:2),]),1,function(u) mean((Ztrunc[1,]-u)^2,na.rm=TRUE))
     p.init[is.na(p.init)]=min(p.init,na.rm=TRUE)
     p.init[p.init==Inf]=min(p.init,na.rm=TRUE)
     p.init=p.init/p.init[round(Ncompet/2)]
     p.ref=p.init[round(Ncompet/2)]
     
+    
     f=function(par){
-      par1=c(par[1:(round(Ncompet/2)-1)],p.ref,par[round(Ncompet/2):(Ncompet-1)])	
+      if(Ncompet==1){ par1=1 } else {
+      	par1=c(par[1:(round(Ncompet/2)-1)],p.ref,par[round(Ncompet/2):(Ncompet-1)])	
+      }
       par1=par1/sum(par1)
       par2=par[Ncompet:length(par)]
-      return(-loglik(Ztrunc[series0,],Ztrunc[seriescompet,],p=par1,sigma=par2,penal.lambda=lambda))
+      return(-loglik(Ztrunc[series0,],rbind(Ztrunc[seriescompet,]),
+      	p=par1,sigma=par2,penal.lambda=lambda))
     }
     par.init=c(rep(1,Ncompet-1),rep(10^2,Ncompet))
     lambda=0
+    
+    
     f(par.init)
-    OPT=optim(par.init, f , method="L-BFGS-B", lower=rep(0,length(par.init)),
+    
+    if(length(par.init)==1) {
+    	OPT=optimize(f,c(0,10^4))
+    	OPT$p=1
+    	OPT$sigma=OPT$minimum
+	    crit=NULL
+	    removeCV=NULL
+    } else {
+    
+ 	   OPT=optim(par.init, f , method="L-BFGS-B", lower=rep(0,length(par.init)),
               upper=rep(Inf,length(par.init)),control=list(factr=10^(-12)))
-    OPT$p=pmax(0,c(OPT$p[1:(round(Ncompet/2)-1)],p.ref,OPT$p[round(Ncompet/2):(Ncompet-1)]))
-    ##OPT$p=pmax(0,c(p.ref,OPT$par[1:(Ncompet-1)]))
-    OPT$p=OPT$p/sum(OPT$p)
-    OPT$sigma=OPT$par[Ncompet:length(OPT$par)]
-    par.init=OPT$par
+              
+              
+   	 OPT$p=pmax(0,c(OPT$p[1:(round(Ncompet/2)-1)],p.ref,OPT$p[round(Ncompet/2):(Ncompet-1)]))
+   	 ##OPT$p=pmax(0,c(p.ref,OPT$par[1:(Ncompet-1)]))
+   	 OPT$p=OPT$p/sum(OPT$p)
+   	 OPT$sigma=OPT$par[Ncompet:length(OPT$par)]
+   	 par.init=OPT$par
     f(par.init)
     crit=NULL
     for(lambda in lambda.series){
@@ -179,7 +202,9 @@ compare_slope_generic_predict_f = function(list_global, model1 = FALSE)
     ##OPT$p=pmax(0,c(p.ref,OPT$par[1:(Ncompet-1)]))
     OPT$p=OPT$p/sum(OPT$p)
     OPT$sigma=OPT$par[Ncompet:length(OPT$par)]
-    OPT
+    
+    }
+    
   }