EBSeq/R/LogNMultiEMP.R

   1 LogNMulti <-
   2 function(Input, InputSP, EmpiricalR, EmpiricalRSP, NumOfEachGroup, AlphaIn, BetaIn,  PIn, NoneZeroLength, AllParti, Conditions)
   3 {
   4
   5         #For each gene (m rows of Input---m genes)
   6         #Save each gene's F0, F1 for further likelihood calculation.
   7                 FList=sapply(1:nrow(AllParti),function(i)sapply(1:nlevels(as.factor(AllParti[i,])),
   8                                    function(j)f0(do.call(cbind,InputSP[AllParti[i,]==j]),AlphaIn, BetaIn,
   9                                                                          do.call(cbind,EmpiricalRSP[AllParti[i,]==j]), NumOfEachGroup, log=T)),
  10                                           simplify=F)
  11                 FPartiLog=sapply(FList,rowSums)
  12                 FMat=exp(FPartiLog)
  13                 rownames(FMat)=rownames(Input)
  14         #Get z
  15                 #Use data.list in logfunction
  16         PInMat=matrix(rep(1,nrow(Input)),ncol=1)%*%matrix(PIn,nrow=1)
  17                 FmultiP=FMat*PInMat
  18                 Denom=rowSums(FmultiP)
  19                 ZEach=apply(FmultiP,2,function(i)i/Denom)
  20                 zNaNName1=names(Denom)[is.na(Denom)]
  21                 # other NAs in LikeFun
  22                 LF=ZEach*(log(FmultiP))
  23                 zNaNMore=rownames(LF)[which(is.na(rowSums(LF)))]
  24                 zNaNName=unique(c(zNaNName1,zNaNMore))
  25                 zGood=which(!rownames(LF)%in%zNaNName)
  26                 ZEachGood=ZEach[zGood,]
  27                 ###Update P
  28         PFromZ=colSums(ZEach[zGood,])/length(zGood)
  29         FGood=FMat[zGood,]
  30                 ### MLE Part ####
  31         # Since we dont wanna update p and Z in this step
  32         # Each Ng for one row
  33
  34                 NumGroupVector=rep(c(1:NoneZeroLength),NumOfEachGroup)
  35
  36                 NumGroupVector.zGood=NumGroupVector[zGood]
  37                 NumOfEachGroup.zGood=tapply(NumGroupVector.zGood,NumGroupVector.zGood,length)
  38
  39         StartValue=c(AlphaIn, BetaIn,PIn[-1])
  40                 InputSPGood=sapply(1:length(InputSP),function(i)InputSP[[i]][zGood,],simplify=F)
  41         EmpiricalRSPGood=sapply(1:length(EmpiricalRSP),function(i)EmpiricalRSP[[i]][zGood,],simplify=F)
  42
  43                 Result<-optim(StartValue,LikefunMulti,InputPool=list(InputSPGood,Input[zGood,],ZEach[zGood,],
  44                                          NoneZeroLength,EmpiricalR[zGood, ],EmpiricalRSPGood, NumOfEachGroup.zGood, AllParti))
  45                 AlphaNew= Result$par[1]
  46                 BetaNew=Result$par[2:(1+NoneZeroLength)]
  47         PNewNo1=Result$par[(2+NoneZeroLength):length(Result$par)]
  48                 PNew=c(1-sum(PNewNo1),PNewNo1)
  49                 ##
  50         Output=list(AlphaNew=AlphaNew,BetaNew=BetaNew,PNew=PNew,ZEachNew=ZEach, ZEachGood=ZEachGood,
  51                                         PFromZ=PFromZ, zGood=zGood, zNaNName=zNaNName,FGood=FGood)
  52         Output
  53     }
  54