EBSeq/R/LogN.R

   1 LogN <-
   2 function(Input, InputSP, EmpiricalR, EmpiricalRSP, NumOfEachGroup, AlphaIn, BetaIn,  PIn, NoneZeroLength)
   3 {
   4     #2 condition case (skip the loop then maybe run faster? Code multi condition cases later)
   5
   6         #For each gene (m rows of Input---m genes)
   7         #Save each gene's F0, F1 for further likelihood calculation.
   8
   9         #Get F0 for EE
  10         F0=f0(Input,  AlphaIn, BetaIn, EmpiricalR, NumOfEachGroup, log=F)
  11         #Get F1 for DE
  12         F1=f1(InputSP[[1]], InputSP[[2]], AlphaIn, BetaIn, EmpiricalRSP[[1]],EmpiricalRSP[[2]], NumOfEachGroup, log=F)
  13
  14         #Get z
  15                 #Use data.list in logfunction
  16
  17                 z.list=PIn*F1/(PIn*F1+(1-PIn)*F0)
  18                 zNaNName=names(z.list)[is.na(z.list)]
  19                 zGood=which(!is.na(z.list))
  20                 ###Update P
  21         #PFromZ=sapply(1:NoneZeroLength,function(i) sum(z.list[[i]])/length(z.list[[i]]))
  22         PFromZ=sum(z.list[zGood])/length(z.list[zGood])
  23         F0Good=F0[zGood]
  24                 F1Good=F1[zGood]
  25                 ### MLE Part ####
  26         # Since we dont wanna update p and Z in this step
  27         # Each Ng for one row
  28
  29                 NumGroupVector=rep(c(1:NoneZeroLength),NumOfEachGroup)
  30
  31                 NumGroupVector.zGood=NumGroupVector[zGood]
  32                 NumOfEachGroup.zGood=tapply(NumGroupVector.zGood,NumGroupVector.zGood,length)
  33
  34         StartValue=c(AlphaIn, BetaIn,PIn)
  35
  36                 Result<-optim(StartValue,Likefun,InputPool=list(InputSP[[1]][zGood,],InputSP[[2]][zGood,],Input[zGood,],z.list[zGood], NoneZeroLength,EmpiricalR[zGood, ],EmpiricalRSP[[1]][zGood,], EmpiricalRSP[[2]][zGood,], NumOfEachGroup.zGood))
  37         #LikeOutput=Likelihood( StartValue, Input , InputSP , PNEW.list, z.list)
  38                 AlphaNew= Result$par[1]
  39                 BetaNew=Result$par[2:(1+NoneZeroLength)]
  40         PNew=Result$par[2+NoneZeroLength]
  41                 ##
  42         Output=list(AlphaNew=AlphaNew,BetaNew=BetaNew,PNew=PNew,ZNew.list=z.list,PFromZ=PFromZ, zGood=zGood, zNaNName=zNaNName,F0Out=F0Good, F1Out=F1Good)
  43         Output
  44     }
  45