##############################################################################################

D.make=function(nQTL,type="RI",a=1:nQTL,d=1:nQTL,aa=NULL,dd=NULL,ad=NULL){
  if(type=="BC"){
    D.matrix=matrix(c(0.5,-0.5),2,1)
    row.names(D.matrix)=c(2,1)
    colnames(D.matrix)="a1"
    if(nQTL>1){
      for(i in 2:nQTL){
        D.matrix=rbind(cbind(0.5,D.matrix),cbind(-0.5,D.matrix))
      }
      colnames(D.matrix)=paste("a",1:nQTL,sep = "")
      D2=matrix(0,2^nQTL,nQTL)
      for(i in 1:nQTL){
        D2[,i]=rep(rep(c(2,1),each=2^(nQTL-i)),2^(i-1))
      }
      row.names(D.matrix)=apply(D2,1,paste,collapse="")
    }
  } else {
    D.matrix=matrix(c(1,0,-1,-0.5,0.5,-0.5),3,2)
    row.names(D.matrix)=c(2,1,0)
    colnames(D.matrix)=c("a1","d1")
    if(nQTL>1){
      for(i in 2:nQTL){
        D.matrix=rbind(cbind(1,-0.5,D.matrix),cbind(0,0.5,D.matrix),cbind(-1,-0.5,D.matrix))
      }
      colnames(D.matrix)=paste(rep(c("a","d"),nQTL),rep(1:nQTL,each=2),sep = "")
      D2=matrix(0,3^nQTL,nQTL)
      for(i in 1:nQTL){
        D2[,i]=rep(rep(c(2,1,0),each=3^(nQTL-i)),3^(i-1))
      }
      row.names(D.matrix)=apply(D2,1,paste,collapse="")
    }
  }

  if(nQTL>1){
    if(!is.null(aa)){
      if(class(aa)=="character"){
        aa=t(combn(nQTL,2))
      } else {aa=matrix(aa,length(aa)/2,2)}
      for(j in 1:nrow(aa)){
        D1=matrix(D.matrix[,colnames(D.matrix)==paste("a",aa[j,1],sep = "")]*
                    D.matrix[,colnames(D.matrix)==paste("a",aa[j,2],sep = "")],nrow(D.matrix),1)
        colnames(D1)=paste("a",aa[j,1],":","a",aa[j,2],sep = "")
        D.matrix=cbind(D.matrix,D1)
      }
    }
    if(type!="BC"){
      if(!is.null(dd)){
        if(class(dd)=="character"){
          dd=t(combn(nQTL,2))
        } else {dd=matrix(dd,length(dd)/2,2)}
        for(j in 1:nrow(dd)){
          D1=matrix(D.matrix[,colnames(D.matrix)==paste("d",dd[j,1],sep = "")]*
                      D.matrix[,colnames(D.matrix)==paste("d",dd[j,2],sep = "")],nrow(D.matrix),1)
          colnames(D1)=paste("d",dd[j,1],":","d",dd[j,2],sep = "")
          D.matrix=cbind(D.matrix,D1)
        }
      }
      if(!is.null(ad)){
        if(class(ad)=="character"){
          ad=t(combn(nQTL,2))
          ad=rbind(ad,cbind(ad[,2],ad[,1]))
          ad=ad[order(ad[,1],ad[,2]),]
        } else {ad=matrix(ad,2,length(ad)/2)}
        for(j in 1:nrow(ad)){
          D1=matrix(D.matrix[,colnames(D.matrix)==paste("a",ad[j,1],sep = "")]*
                      D.matrix[,colnames(D.matrix)==paste("d",ad[j,2],sep = "")],nrow(D.matrix),1)
          colnames(D1)=paste("a",ad[j,1],":","d",ad[j,2],sep = "")
          D.matrix=cbind(D.matrix,D1)
        }
      }
    }
    a=a[a<=nQTL]
    a=unique(a)
    a=sort(a)
    a.cut=(1:nQTL)[!(1:nQTL)%in%a]

    d=d[d<=nQTL]
    d=unique(d)
    d=sort(d)
    d.cut=(1:nQTL)[!(1:nQTL)%in%d]
    
    DD=D.matrix
    
    if(length(a.cut)>0 & length(d.cut)>0){
      DD=D.matrix[,-c((a.cut*2-1),(d.cut*2))]
    } else {
      if(length(a.cut)>0){
        d.cut=NA
        DD=D.matrix[,-(a.cut*2-1)]
      }
      if(length(d.cut)>0){
        a.cut=NA
        DD=D.matrix[,-(d.cut*2)]
      }
    }
    
    if(class(DD)=="numeric"){
      cut0=c((a.cut*2-1),(d.cut*2))
      cut0=cut0[!is.na(cut0)]
      DD=matrix(DD,nrow(D.matrix),1)
      colnames(DD)=colnames(D.matrix)[setdiff((1:ncol(D.matrix)), cut0)]
      row.names(DD)=row.names(D.matrix)
    } else if (ncol(DD)==0){
      DD=1
    }
    
   D.matrix=DD

  } else {
    if(sum(d,na.rm = T)==0){
      DD=matrix(D.matrix[,1],nrow(D.matrix),1)
      colnames(DD)="a1"
      row.names(DD)=row.names(D.matrix)
      D.matrix=DD
    }
    if(sum(a,na.rm = T)==0){
      DD=matrix(D.matrix[,2],nrow(D.matrix),1)
      colnames(DD)="d1"
      row.names(DD)=row.names(D.matrix)
      D.matrix=DD
    }
    if(sum(c(a,d),na.rm = T)==0){
      D.matrix=1
    }
  }

  return(D.matrix)
}


################################################################################################
Q.make=function(QTL,marker,geno=NULL,interval=F,cM=T,type="RI",ng=2){
  # marker: n.marker*2 matrix
  # geno: n.ind*n.marker matrix

  cr=QTL[,1]
  QTL=QTL[,2]
  
  if(cM==T){
    QTL=QTL/100
    marker[,2]=marker[,2]/100
  }
  
  nd=3
  if(type=="BC"){
    nd=2
    type.fun=function(d1,d2,ng){
      MN=bcp2M(d1+d2,ng)
      A=bcp3M(d1,d2,ng)
      Q0=c(A[1]/MN[1],A[3]/MN[1],
           A[2]/MN[2],A[7]/MN[2],
           A[4]/MN[3],A[5]/MN[3],
           A[6]/MN[4],A[8]/MN[4])
      Q1=matrix(Q0,4,2,byrow = T)
      colnames(Q1)=c("QQ","Qq")
      rownames(Q1)=c(22,21,12,11)
      return(Q1)
    }
  } else if (type=="RI"){
    type.fun=function(d1,d2,ng){
      MN=RI2(d1+d2,ng)
      A=RI3(d1,d2,ng)
      Q0=c(A[1:3]/MN[1],c(A[4],A[5]+A[6],A[7])/MN[2],A[8:10]/MN[3],
           c(A[11],A[12]+A[13],A[14])/MN[4],
           c(A[15]+A[16],sum(A[17:20]),A[15]+A[16])/MN[5],
           c(A[14],A[12]+A[13],A[11])/MN[6],
           A[10:8]/MN[7],c(A[7],A[5]+A[6],A[4])/MN[8],A[3:1]/MN[9])
      Q1=matrix(Q0,9,3,byrow = T)
      colnames(Q1)=c("QQ","Qq","qq")
      rownames(Q1)=c(22,21,20,12,11,10,"02","01","00")
      return(Q1)
    }
  } else if (type=="AI"){
    type.fun=function(d1,d2,ng){
      MN=AI2(d1+d2,ng)
      A=AI3(d1,d2,ng)
      Q0=c(A[1:3]/MN[1],c(A[4],A[5]+A[6],A[7])/MN[2],A[8:10]/MN[3],
           c(A[11],A[12]+A[13],A[14])/MN[4],
           c(A[15]+A[16],sum(A[17:20]),A[15]+A[16])/MN[5],
           c(A[14],A[12]+A[13],A[11])/MN[6],
           A[10:8]/MN[7],c(A[7],A[5]+A[6],A[4])/MN[8],A[3:1]/MN[9])
      Q1=matrix(Q0,9,3,byrow = T)
      colnames(Q1)=c("QQ","Qq","qq")
      rownames(Q1)=c(22,21,20,12,11,10,"02","01","00")
      return(Q1)
    }
  }
  
  marker=cbind(marker,1:nrow(marker))
  
  cr0=marker[marker[,1]==cr[1],]
  
  marker1=cr0[max(which(cr0[,2]<QTL[1])),]
  marker2=cr0[max(which(cr0[,2]<QTL[1]))+1,]
  if(interval==T){
    marker2=cr0[min(which(cr0[,2]>QTL[1])),]
  }
  marker0=as.numeric(c(marker1[3],marker2[3]))
  
  d1=as.numeric(QTL[1]-marker1[2])
  d2=as.numeric(marker2[2]-QTL[1])
  Q1=type.fun(d1,d2,ng)
  
  Q1=list(Q1)
  m=1
  
  n=length(QTL)
  if(n>1){
    for(m in 2:n){
      cr0=marker[marker[,1]==cr[m],]
      
      marker1=cr0[max(which(cr0[,2]<QTL[m])),]
      marker2=cr0[max(which(cr0[,2]<QTL[m]))+1,]
      if(interval==T){
        marker2=cr0[min(which(cr0[,2]>QTL[m])),]
      }
      marker0=as.numeric(c(marker0,marker1[3],marker2[3]))
      
      d1=as.numeric(QTL[m]-marker1[2])
      d2=as.numeric(marker2[2]-QTL[m])
      Q2=type.fun(d1,d2,ng)
      Q1[[m]]=Q2
    }
  }
  

  if(!is.null(geno)){
    red.genotype=geno[,marker0]
    cp.matrix=matrix(0,nrow(geno),nd^n)
    D2=matrix(0,nd^n,n)
    for(i in 1:n){
      D2[,i]=rep(rep(2:(3-nd),each=nd^(n-i)),nd^(i-1))
    }
    qname=apply(D2,1,paste,collapse="")
    for(j in 1:nrow(geno)){
      geno.j=red.genotype[j,]
      pq=c()
      for(k in 1:length(QTL)){
        g0=geno.j[(k*2-1):(k*2)]
        if(is.na(g0[1]) & is.na(g0[2])){
          a=Q1[[k]]
          a=matrix(unlist(a),nrow(a),ncol(a))
          pq0=apply(a,2,mean)
        } else if (!is.na(g0[1]) & is.na(g0[2])){
          a=row.names(Q1[[k]])
          a=as.numeric(unlist(strsplit(a,split="",fixed=T)) )
          a=t(matrix(a,2))
          a=Q1[[k]][a[,1]==as.numeric(g0[1]),]
          a=matrix(unlist(a),nrow(a),ncol(a))
          pq0=apply(a,2,mean)
        } else if (is.na(g0[1]) & !is.na(g0[2])){
          a=row.names(Q1[[k]])
          a=as.numeric(unlist(strsplit(a,split="",fixed=T)) )
          a=t(matrix(a,2))
          a=Q1[[k]][a[,2]==as.numeric(g0[2]),]
          a=matrix(unlist(a),nrow(a),ncol(a))
          pq0=apply(a,2,mean)
        } else {
          g0=paste(g0[1],g0[2],sep = "")
          pq0=Q1[[k]]
          pq0=pq0[row.names(pq0)==g0]
        }
        pq=rbind(pq,pq0)
      }
      
      for(i in 1:nd^n){
        a=qname[i]
        a=-as.numeric(unlist(strsplit(a,split="",fixed=T)))+3
        pq1=pq[1:length(QTL),a]
        if(is.matrix(pq1)){cp.matrix[j,i]=prod(as.numeric(diag(pq1)))
        } else {cp.matrix[j,i]=as.numeric(pq1)}
      }
    }
    colnames(cp.matrix)=qname
    Q1[[m+1]]=cp.matrix
  }
  
  names(Q1)[1:m]=c(paste("Q.matrix",1:m,sep = ""))
  if(length(Q1)==m+1){names(Q1)[m+1]="cp.matrix"}
  return(Q1)
}


################################## MIM ##############################################

EM.MIM=function(D.matrix,cp.matrix,y,E.vector=NULL,X=NULL,beta=NULL,variance=NULL,conv=10^-5,console=F){
  Y=y[!is.na(y)]
  cp.matrix=cp.matrix[!is.na(y),]
  ind=length(Y)
  g=nrow(D.matrix)
  eff=ncol(D.matrix)
  
  if(is.null(E.vector)){E.vector=rep(0,eff)}
  if(is.null(X)){
    X=matrix(1,ind,1)
    if(length(beta>1)){
      beta=mean(beta)
    }
  }else if(is.vector(X)){
    X=matrix(X,length(X),1)
    X=X[!is.na(y),]
  }
  if(is.null(beta)){
    beta=matrix(rep(mean(Y),ncol(X)),ncol(X),1)
  }else if(is.numeric(beta)){
    beta=matrix(rep(beta,ncol(X)),ncol(X),1)
  }
  if(is.null(variance)){variance=var(Y)}
  sigma=sqrt(variance)
  
  Delta=1
  time=0
  
  while (max(abs(Delta))>conv) {
    muji.matrix=t(D.matrix%*%E.vector%*%matrix(1,1,ind))+X%*%beta%*%matrix(1,1,g)
    
    PI.matrix=matrix(0,ind,g)
    for(j in 1:ind){
      P0=c()
      for(i in 1:g){
        P0[i]=cp.matrix[j,i]*dnorm((Y[j]-muji.matrix[j,i])/sigma,0,1)
      }
      if(sum(P0)!=0){P0=P0/sum(P0)
      } else {P0=rep(1/g,g)}
      PI.matrix[j,]=P0
    }
    
    r.vector=c()
    for(i in 1:eff){
      r01=t(Y-X%*%beta)%*%PI.matrix%*%D.matrix[,i]
      r02=matrix(1,1,ind)%*%PI.matrix%*%(D.matrix[,i]*D.matrix[,i])
      r.vector[i]=r01/r02
    }
    
    M.matrix=matrix(0,eff,eff)
    for(i in 1:eff){
      for(j in 1:eff){
        if(i!=j){
          M01=matrix(1,1,ind)%*%PI.matrix%*%(D.matrix[,i]*D.matrix[,j])
          M02=matrix(1,1,ind)%*%PI.matrix%*%(D.matrix[,i]*D.matrix[,i])
          M.matrix[i,j]=M01/M02
        }
      }
    }
    
    E.t=r.vector-M.matrix%*%E.vector
    
    beta.t=solve(t(X)%*%X)%*%t(X)%*%(Y-PI.matrix%*%D.matrix%*%E.t)
    
    V.matrix=matrix(0,eff,eff)
    for(i in 1:eff){
      for(j in 1:eff){
        V.matrix[i,j]=matrix(1,1,ind)%*%PI.matrix%*%(D.matrix[,i]*D.matrix[,j])
      }
    }
    
    
    sigma.t=sqrt((t(Y-X%*%beta.t)%*%(Y-X%*%beta.t)-t(Y-X%*%beta.t)%*%PI.matrix%*%D.matrix%*%E.t*2+t(E.t)%*%V.matrix%*%E.t)/ind)
    
    Delta=E.t-E.vector
    time=time+1
    if(console==T){
      Ep=round(E.t,5)
      sp=round(sigma.t^2,5)
      print(c("Time",time,"effect",Ep,"var",sp))
    }
    
    E.vector=E.t
    beta=beta.t
    sigma=sigma.t
  }
  
  PI.matrix=matrix(0,ind,g)
  for(j in 1:ind){
    P0=c()
    for(i in 1:g){
      P0[i]=cp.matrix[j,i]*dnorm((Y[j]-muji.matrix[j,i])/sigma,0,1)
    }
    P0=P0/sum(P0)
    PI.matrix[j,]=P0
  }
  E.vector=c(E.vector)
  names(E.vector)=colnames(D.matrix)
  colnames(PI.matrix)=colnames(cp.matrix)
  variance=sigma^2
  
  L0=c()
  L1=c()
  for(k in 1:nrow(cp.matrix)){
    L00=c()
    L01=c()
    for(m in 1:nrow(D.matrix)){
      L00[m]=cp.matrix[k,m]*dnorm((Y[k]-mean(X%*%beta))/sigma)
      L01[m]=cp.matrix[k,m]*dnorm((Y[k]-(mean(X%*%beta)+D.matrix[m,]%*%E.vector))/sigma)
    }
    L0[k]=sum(L00)
    L1[k]=sum(L01)
  }
  LRT=-2*sum(log(L0[!is.na(L0) & !is.na(L1)]/L1[!is.na(L0) & !is.na(L1)]))
  
  result=list(E.vector=E.vector,beta=as.numeric(beta),variance=as.numeric(variance),PI.matrix=PI.matrix,LRT=LRT,iteration.time=time)
  return(result)
}


EM.MIM2=function(QTL,marker,geno,D.matrix,cp.matrix=NULL,y,yu=NULL,sele.g="n",tL=NULL,tR=NULL,cM=T,type="RI",ng=2,
                 E.vector=NULL,X=NULL,beta=NULL,variance=NULL,conv=10^-5,console=F){
  Y=y[!is.na(y)]
  Yu=yu[!is.na(yu)]
  if(sele.g=="p" | sele.g=="pf"){
    Ya=c(y,yu)
  }else{Ya=y}
  geno=geno[which(!is.na(y)),]
  if(!is.null(cp.matrix)){cp.matrix=cp.matrix[which(!is.na(Ya)),]}
  if(is.null(E.vector)){E.vector=rep(0,ncol(D.matrix))}
  if(is.null(X)){
    X=matrix(1,length(Ya),1)
    if(length(beta>1)){
      beta=mean(beta)
    }
  }else if(is.vector(X)){
    X=matrix(X,length(Ya),1)
  }
  X=as.matrix(X[which(!is.na(Ya)),])
  if(is.null(beta)){
    beta=matrix(rep(mean(Y),ncol(X)),ncol(X),1)
  }else if(is.numeric(beta)){
    beta=matrix(rep(beta,ncol(X)),ncol(X),1)
  }
  nQTL=nrow(QTL)
  
  if(sele.g=="p"){
    if(nQTL>1){QTL=QTL[order(QTL[,1],QTL[,2]),]}
    marker=marker[order(marker[,1],marker[,2]),]
    result=QTL.p(QTL,marker,geno,D.matrix, ys=Y, yu=Yu, cM=cM, type=type, ng=ng, 
                 E.vector=E.vector, X=X, beta=beta, variance=variance, 
                 cp.matrix=cp.matrix, conv=conv, console=console)
    model="proposed model of selective genotyping"
    result[[8]]=model
    names(result)[8]="model"
  } else if (sele.g=="pf"){
    if(nQTL>1){QTL=QTL[order(QTL[,1],QTL[,2]),]}
    marker=marker[order(marker[,1],marker[,2]),]
    result=QTL.p(QTL,marker,geno,D.matrix, ys=Y, yu=Yu, cM=cM, type=type, ng=ng, 
                 E.vector=E.vector, X=X, beta=beta, variance=variance, 
                 cp.matrix=cp.matrix, conv=conv, console=console, pf = T)
    model="population frequency-based model of selective genotyping"
    result[[8]]=model
    names(result)[8]="model"
  }else if(sele.g=="t"){
    if(nQTL>1){QTL=QTL[order(QTL[,1],QTL[,2]),]}
    marker=marker[order(marker[,1],marker[,2]),]
    if(is.null(tL)){tL=min(Yu)}
    if(is.null(tR)){tR=max(Yu)}
    result=QTL.t(QTL, marker, geno, D.matrix, ys=Y, yu=yu, tL, tR, cM=cM, type=type, ng=ng, 
                 E.vector=E.vector, beta=beta,  X=X, variance=variance, 
                 cp.matrix=cp.matrix, conv=conv, console=console)
    model="truncated model of selective genotyping"
    result[[8]]=model
    names(result)[8]="model"
  }else{
    if(is.null(cp.matrix)){
      cp.matrix=Q.make(QTL,marker,geno,cM=cM,type=type,ng=ng)[[(nrow(QTL)+1)]]
    }
    re=EM.MIM(D.matrix,cp.matrix,Y,E.vector=E.vector,X=X,beta=beta,variance=variance,conv=conv,console=console)
    result=list(QTL=QTL,E.vector=re[[1]],beta=re[[2]],variance=re[[3]],PI.matrix=re[[4]],LRT=re[[5]],iteration.time=re[[6]],model="complete genotyping model")
  }
  names(result[[2]])=colnames(D.matrix)
  return(result)
}


############################### IM #################################

detectQTL=function(effect,LRT.threshold,minQTLdist=10){
  LRT=effect[,c(1,2,ncol(effect))]
  LRT[LRT[,3]<LRT.threshold,3]=0
  det0=c()
  for(i in unique(LRT[,1])){
    LRTcr=LRT[LRT[,1]==i,2:3]
    detcr=rep(0,nrow(LRTcr))
    if(sum(LRTcr[,2])>0){
      for(j in 1:nrow(LRTcr)){
        LRTdet=LRTcr[LRTcr[,1]>=max(c(min(LRTcr[,1]),LRTcr[j,1]-minQTLdist)) & 
                       LRTcr[,1]<=min(c(max(LRTcr[,1]),LRTcr[j,1]+minQTLdist)),2]
        LRT0=LRTcr[j,2]
        if(LRT0==max(LRTdet) & LRT0>0){detcr[j]=1
        }else{detcr[j]=0}
      }
    }
    det0=c(det0,detcr)
  }
  detect.QTL=effect[det0==1,]
  detect.QTL
}


EM.IM=function(marker,geno,y,method="EM",type="RI",ng=2,D.matrix=NULL,cM=T,speed=1,conv = 10^-5,d.eff=T,
               LRT.thre=T,simu=1000,lv=0.05,detect=T,minQTLdist=15,chart = "chr",console = F){
  Y=y
  
  if(cM==F){
    marker[,2]=marker[,2]*100
  }
  
  if(is.null(D.matrix)){
    if(type=="BC"){
      D.matrix=matrix(c(0.5,-0.5),2,1)
      row.names(D.matrix)=c(2,1)
      colnames(D.matrix)="a1"
    } else {
      a1=c(1,0,-1)
      d1=c(-0.5,0.5,-0.5)
      if(d.eff==T){
        D.matrix=cbind(a1,d1)
        row.names(D.matrix)=c(2,1,0)
      } else {
        D.matrix=matrix(a1,3,1)
        row.names(D.matrix)=c(2,1,0)
        colnames(D.matrix)="a1"
      }
    }
  }

  if(method=="EM"){
    meth=function(D.matrix,cp.matrix,Y,conv){
      EM=EM.MIM(D.matrix,cp.matrix,Y,conv = conv)
      eff=as.numeric(EM$E.vector)
      mu0=as.numeric(EM$beta)
      sigma=sqrt(as.numeric(EM$variance))
      result=list(eff,mu0,sigma)
      return(result)
    }
  } else if (method=="REG"){
    meth=function(D.matrix,cp.matrix,Y,conv){
      X=cp.matrix%*%D.matrix
      fit=lm(Y~X)
      eff=as.numeric(fit$coefficients[-1])
      mu0=as.numeric(fit$coefficients[1])
      ms=anova(fit)$`Mean Sq`
      sigma=ms[2]^0.5
      result=list(eff,mu0,sigma)
      return(result)
    }
  }

  effect=c()
  cr0=unique(marker[,1])
  for(i in cr0){
    cr=marker[marker[,1]==i,]
    for(j in seq(ceiling(floor(min(cr[,2]))+speed),(max(cr[,2])),speed)){
      Q.matrix=Q.make(matrix(c(i,j),1,2),marker,geno,type=type,ng=ng)
      cp.matrix=Q.matrix$cp.matrix
      
      result=meth(D.matrix,cp.matrix,Y,conv)
      eff=result[[1]]
      mu0=result[[2]]
      sigma=result[[3]]

      L0=c()
      L1=c()
      for(k in 1:nrow(cp.matrix)){
        L00=c()
        L01=c()
        for(m in 1:nrow(D.matrix)){
          L00[m]=cp.matrix[k,m]*dnorm((Y[k]-mu0)/sigma)
          L01[m]=cp.matrix[k,m]*dnorm((Y[k]-(mu0+D.matrix[m,]%*%eff))/sigma)
        }
        L0[k]=sum(L00)
        L1[k]=sum(L01)
      }
      
      LRT=-2*sum(log(L0[!is.na(L0) & !is.na(L1)]/L1[!is.na(L0) & !is.na(L1)]))
      
      
      
      eff0=c(i,j,eff,LRT)
      effect=rbind(effect,eff0)
      if(console==T){print(c(paste("chr",i),paste(j,"cM"),paste("LRT",LRT)))}
    }
  }
  row.names(effect)=1:nrow(effect)
  colnames(effect)=c("chr","cM",colnames(D.matrix),"LRT")
  effect=data.frame(effect)
  effect[effect[,ncol(effect)]==Inf,3:ncol(effect)]=0
  
  
  LRT.threshold=NULL
  if(LRT.thre==T){
    LRT.threshold=LRTthre(marker,type,ng,length(Y),25,cM,d.eff,simu,speed,lv,console)
  } else if (class(LRT.thre)=="numeric"){
    LRT.threshold=LRT.thre
  }
  
  
  detect.QTL=NULL
  if(detect==T & class(LRT.threshold)=="numeric"){
    detect.QTL=detectQTL(effect,LRT.threshold,minQTLdist)
  }
  
  if(chart=="chr"){
    for(k in cr0){
      par(mar=c(5,5,4,2))
      par(fig=c(0,1,0.35,1))
      plot(effect[effect$chr==k,]$cM,effect[effect$chr==k,]$LRT,main="",ylab="LRT statistic",
           xlab="",type="l",ylim=c(0,max(c(effect$LRT,LRT.threshold),na.rm = T)),bty="l")
      if(!is.null(LRT.threshold)){abline(h=LRT.threshold,col="red",lwd=2)}
      par(mar=c(2,5,4,2))
      par(fig=c(0,1,0,0.45),new=T)
      plot(effect[effect$chr==k,]$cM,effect[effect$chr==k,]$a1,ylab="effect",main=paste("Chromosome",k),
           xlab="",type="l",ylim=c(min(c(effect$a1,0))*1.4,max(c(effect$a1,0))*1.2),xaxt="n",bty="n")
      abline(h=0,col="black",lwd=2)
    }
  } else if (chart=="all") {
    t=table(effect$chr)
    ncr=cumsum(t)
    par(mar=c(5,5,4,2))
    par(fig=c(0,1,0.3,1))
    plot(effect$LRT,main="",ylab="LRT statistic",xlab="",type="n",
         ylim=c(0,max(c(effect$LRT,LRT.threshold),na.rm = T)),xaxt="n",bty="l")
    k0=0
    for(k in cr0){
      cMk=max(effect[effect$chr==k,]$cM)
      points(effect[effect$chr==k,]$cM+k0,effect[effect$chr==k,]$LRT,main="",ylab="LRT statistic",xlab="",
             type="l",ylim=c(0,max(c(effect$LRT,LRT.threshold),na.rm = T)),xaxt="n",bty="l")
      k0=k0+cMk
      abline(v=k0,col="blue",lwd=2)
    }
    abline(v=k0,col="white",lwd=2)
    axis(1,paste("chr",(1:length(ncr))),at=ncr-t/2)
    if(!is.null(LRT.threshold)){abline(h=LRT.threshold,col="red",lwd=2)}
    
    par(mar=c(2,5,4,2))
    par(fig=c(0,1,0,0.5),new=T)
    plot(effect$a1,main="",ylab="effect",xlab="",type="n",
         ylim=c(min(c(effect$a1,0))*1.2,max(c(effect$a1,0))*1.2),xaxt="n",bty="n")
    k0=0
    for(k in cr0){
      cMk=max(effect[effect$chr==k,]$cM)
      points(effect[effect$chr==k,]$cM+k0,effect[effect$chr==k,]$a1,main="",ylab="effect",xlab="",type="l",
             ylim=c(min(c(effect$a1,0))*1.2,max(c(effect$a1,0))*1.2),xaxt="n",bty="n")
      k0=k0+cMk
      abline(v=k0,col="blue",lwd=2)
    }
    abline(v=k0,col="white",lwd=2)
    abline(h=0,col="black",lwd=2)
  }
  par(mar = c(5.1,4.1,4.1,2.1))
  par(fig=c(0,1,0,1))
  result=list(effect=effect,LRT.threshold=LRT.threshold,detect.QTL=detect.QTL)
}



######################### simulaton F2 ##########################

simulate.progeny = function(DNA0,distance){
  
  distance.c = distance[-1]-distance[-length(distance)]
  distance.r = (1-exp(-2*distance.c/100))/2
  distance.r[distance.r<0] = 0.5
  
  index = sample(1:2,1)
  DNA1 = DNA0[,index]
  DNA2 = DNA0[,-index]
  
  DNA.new = DNA1[1]
  for(i in 2:nrow(DNA0)){
    x = runif(1)
    if(x < distance.r[i-1]){
      DNA.new[i] = DNA2[i]
      D = DNA1
      DNA1 = DNA2
      DNA2 = D
    } else { DNA.new[i] = DNA1[i] }
  }
  
  return(DNA.new)
}

progeny=function(QTL,marker,type="RI",ng=2,E.vector=NULL,h2=0.5,size=200){
  n=nrow(QTL)
  QTLs=cbind(QTL,0)
  marker=cbind(marker,1:nrow(marker))
  marker0=rbind(as.matrix(marker),as.matrix(QTLs))
  marker0=marker0[order(marker0[,1],marker0[,2]),]
  dis=marker0[,2]
  p1=rep(2,nrow(marker0))
  p2=rep(0,nrow(marker0))
  F1=matrix(c(p1,p2),nrow(marker0),2)
  
  D.matrix=D.make(n,type=type,aa="all",dd="all",ad="all")
  if(is.null(E.vector)){E.vector=rep(1,ncol(D.matrix))
  } else if (length(E.vector)!=ncol(D.matrix)){
    E0=rep(0,ncol(D.matrix))
    for(i in 1:length(E.vector)){
      E0[colnames(D.matrix)%in%(names(E.vector)[i])]=E.vector[i]
      if(!((names(E.vector)[i])%in%colnames(D.matrix))){
        name0=strsplit(names(E.vector)[i],split="")
        name1=c()
        for(k in c(4,5,3,1,2)){
          name1=paste(name1,name0[[1]][k],sep = "")
        }
        names(E.vector)[i]=name1
        E0[colnames(D.matrix)%in%(names(E.vector)[i])]=E.vector[i]
      }
    }
    E.vector=E0
  }
  names(E.vector)=colnames(D.matrix)
  
  
  prog.g=list()
  if (type=="BC"){
    p0=cbind(p1,p1)
    for(i in 1:size){
      pr1=simulate.progeny(F1,dis)
      pr2=simulate.progeny(p0,dis)
      prog.g[[i]]=cbind(pr1,pr2)
    }
    if(ng>1){
      for(g in 2:ng){
        prog.g2=list()
        for(i in 1:size){
          sele=sample(1:size,1)
          pr1=simulate.progeny(prog.g[[sele]],dis)
          pr2=simulate.progeny(p0,dis)
          prog.g2[[i]]=cbind(pr1,pr2)
        }
        prog.g=prog.g2
      }
    }
    if(n==1){
      VG=sum((E.vector[1:n])^2/2)
    } else {
      VG=VG.BC(QTL,E.vector)
    }
  } else {
    for(i in 1:size){
      pr1=simulate.progeny(F1,dis)
      pr2=simulate.progeny(F1,dis)
      prog.g[[i]]=cbind(pr1,pr2)
    }
    if(ng>2){
      if(type=="RI"){
        for(g in 3:ng){
          prog.g2=list()
          for(i in 1:size){
            sele=sample(1:size,1)
            pr1=simulate.progeny(prog.g[[sele]],dis)
            pr2=simulate.progeny(prog.g[[sele]],dis)
            prog.g2[[i]]=cbind(pr1,pr2)
          }
          prog.g=prog.g2
        }
      } else if (type=="AI"){
        for(g in 3:ng){
          prog.g2=list()
          for(i in 1:size){
            sele=sample(1:size,2)
            pr1=simulate.progeny(prog.g[[sele[1]]],dis)
            pr2=simulate.progeny(prog.g[[sele[2]]],dis)
            prog.g2[[i]]=cbind(pr1,pr2)
          }
          prog.g=prog.g2
        }
      }
    } else if (ng==1){
      for(i in 1:size){
        prog.g[[i]]=F1
      }
    }

    if(n==1){
      VG=sum((E.vector[(1:n)*2-1])^2/2)+sum((E.vector[(1:n)*2])^2/4)
    } else {
      VG=VG.RI(QTL,E.vector)
    }
  }
  prog=matrix(0,size,length(dis))
  for(i in 1:size){
    prog[i,]=apply(prog.g[[i]],1,mean)
  }
  
  VE=VG*(1-h2)/h2
  phe=c()
  for(i in 1:size){
    genoQ=prog[i,marker0[,3]==0]
    genoQ1=genoQ[1]
    if(n>1){
      for(j in 2:n){
        genoQ1=paste(genoQ1,genoQ[j],sep="")
      }
    }
    eff=D.matrix%*%E.vector
    eff=eff[row.names(eff)==genoQ1,]
    phe[i]=eff+rnorm(1,0,sqrt(VE))
  }
  
  marker.prog=prog[,marker0[,3]!=0]
  QTL.prog=prog[,marker0[,3]==0]
  
  output=list(phe=phe,E.vector=E.vector,marker.prog=marker.prog,QTL.prog=QTL.prog)
  return(output)
}


############################### thre ##################################

LRTthre=function(marker,type="RI",ng=2,ns=200,gv=25,cM=T,d.eff=T,simu=1000,speed=1,lv=0.05,console=F){
  if(cM==T){
    marker[,2]=marker[,2]/100
  }
  
  cr=unique(marker[,1])
  ncr=length(cr)
  thre=c()
  if(type=="BC"){
    thre.fun=thre.BC
  } else if (type=="RI"){
    thre.fun=thre.RI
    if(d.eff==F){
      thre.fun=thre.RIa
    }
  } else if (type=="AI"){
    thre.fun=thre.AI
    if(d.eff==F){
      thre.fun=thre.AIa
    }
  }
  
  for(i in 1:ncr){
    d=marker[marker[,1]==cr[i],2]
    d=d[-1]-d[-length(d)]
    thre0=thre.fun(d,ng,simu,ns,gv,console,cr[i],speed)
    thre=cbind(thre,thre0)
  }
  
  max0=function(x){max(x,na.rm = T)}
  permu.max=apply(thre,1,max0)
  
  thre.all=quantile(permu.max,(1-lv),type=3)
  return(thre.all)
}

############################# other ##################################

mixprop=function(QTL, nu, marker, geno, model,cM=T,type="RI",ng=2,cp.matrix=NULL){ 
  nQTL=nrow(QTL)
  if(nQTL>1){QTL=QTL[order(QTL[,1],QTL[,2]),]}
  marker=marker[order(marker[,1],marker[,2]),]
  
  Q3=3^nQTL
  Q2=2^nQTL
  if(is.null(cp.matrix)){
    cp.matrix=Q.make(QTL,marker,geno,cM=cM,type=type,ng=ng)[[(nQTL+1)]]
  }
  QTL.freq=cp.matrix

  if(cM==T){
    QTL[,2]=QTL[,2]/100
    marker[,2]=marker[,2]/100
  }
  
  if(model==2){
    mix.prop=QTL.freq
    popu.freq=NULL
    freq.u=NULL
  } else {
    N=nrow(QTL.freq)+nu
    freq.s=colSums(QTL.freq)/N
    
    if(nQTL==1){
      if(type=="BC"){
        popu.freq=c(1-0.5^ng,0.5^ng)
        Qn=2
      } else if (type=="RI"){
        hetro=0.5^(ng-1)
        popu.freq=c((1-hetro)/2,hetro,(1-hetro)/2)
        Qn=3
      } else {
        popu.freq=c(0.25,0.5,0.25)
        Qn=3
      }
      
    } else {
      M=c(1,0)
      GAM=matrix(0,Q2,nQTL)
      Pg=c(2,1,0)
      PG=matrix(0,Q3,nQTL)
      rmn=rep(0,nQTL-1)
      for(n0 in 1:nQTL){
        GAM[,n0]=rep(rep(M,2^(n0-1)),each=2^(nQTL-n0))
        PG[,n0]=rep(rep(Pg,3^(n0-1)),each=3^(nQTL-n0))
        if(n0==nQTL){break}
        if(QTL[n0,1]==QTL[(n0+1),1]){
          rmd=QTL[(n0+1),2]-QTL[n0,2]
          rmn[n0]=(1-exp(-2*rmd))/2
        }else{
          rmn[n0]=0.5
        }
      }
      
      if(type=="AI" & ng>2){
        k0=seq(0,(ng-1),2)
        rmn0=matrix(0,(nQTL-1),length(k0))
        for(i in 1:length(k0)){
          rmn0[,i]=choose((ng-1),k0[i])*rmn^k0[i]*(1-rmn)^(ng-1-k0[i])
        }
        rmn=1-apply(rmn0,1,sum)
      }
      
      gamf=GAM.f=matrix(0,Q2,nQTL-1)
      GAM.freq=rep(1,Q2)
      for(nr in 1:(nQTL-1)){
        gamf[,nr]=GAM[,nr]==GAM[,nr+1]
        for(ngam in 1:Q2){
          GAM.f[ngam,nr]=ifelse(gamf[ngam,nr]==1,1-rmn[nr],rmn[nr])
        }
        GAM.freq=GAM.freq*GAM.f[,nr]
      }
      
      if(type=="BC"){
        popu.freq=GAM.freq/2
        if(ng>1){
          G.matrix=matrix(0,Q2,Q2)
          G.matrix[1,1]=1
          G.matrix[Q2,]=GAM.freq/2
          for(i in 2:(Q2-1)){
            p0=apply(matrix(t(GAM[,which(GAM[i,]!=0)])==GAM[i,which(GAM[i,]!=0)],nrow=Q2,byrow=T),1,sum)==sum(GAM[i,]!=0)
            p1=popu.freq[p0]/sum(popu.freq[p0])
            G.matrix[i,p0]=p1
          }
          for(i in 1:(ng-1)){
            popu.freq=crossprod(G.matrix,popu.freq)
          }
        }
        Qn=Q2
      }else{
        
        Gfreq=rep(0,(Q2)*(Q2))
        Gtype=matrix(0,(Q2)*(Q2),nQTL)
        wg=0
        for(ga1 in 1:Q2){
          for(ga2 in 1:Q2){
            wg=wg+1
            Gfreq[wg]=(GAM.freq[ga1]/2)*(GAM.freq[ga2]/2)
            Gtype[wg,]=GAM[ga1,]+GAM[ga2,]
          }
        }
        popu.freq=rep(0,Q3)
        ind=rep(0,(Q2)*(Q2))
        for(pfi in 1:Q3){
          for(gti in 1:((Q2)*(Q2))){
            if(sum(Gtype[gti,]==PG[pfi,])==nQTL){ind[gti]=pfi}
          }
          popu.freq[pfi]=sum(Gfreq[ind==pfi])
        }
        
        if(type=="RI" & ng>2){
          G.matrix=matrix(0,Q3,Q3)
          G.matrix[1,1]=1
          G.matrix[Q3,Q3]=1
          for(i in 2:(Q3-1)){
            p0=apply(matrix(t(PG[,which(PG[i,]!=1)])==PG[i,which(PG[i,]!=1)],nrow=Q3,byrow=T),1,sum)==sum(PG[i,]!=1)
            p1=popu.freq[p0]/sum(popu.freq[p0])
            G.matrix[i,p0]=p1
          }
          for(i in 1:(ng-1)){
            popu.freq=crossprod(G.matrix,popu.freq)
          }
        }
        Qn=Q3
      }
    }
    
    freq.u=as.vector(popu.freq)-freq.s
    if(model==3){
      freq.u=as.vector(popu.freq)
      names(freq.u)=names(freq.s)
    }
    freq.u[freq.u<0]=0
    Freq.u=matrix(rep(freq.u/sum(freq.u),each=nu),nu,Qn)
    mix.prop=rbind(QTL.freq,Freq.u)
  }
  re=list(popu.freq=popu.freq,mix.prop=mix.prop,QTL=QTL,marker=marker,freq.u=freq.u)
  return(re)
}


QTL.p=function(QTL, marker, geno, D.matrix, ys, yu, cM=T, type="RI", ng=2, 
               E.vector=NULL, X=NULL, beta=NULL, variance=NULL, cp.matrix=NULL, conv = 10^-5, console=F, pf=F){
  y=c(ys,yu)
  Y=y[!is.na(y)]
  N=length(Y)
  nu=length(yu[!is.na(yu)])
  nQTL=nrow(QTL)
  model=1
  if(pf==T){model=3}
  mp=mixprop(QTL, nu, marker, geno, model=model, cM=cM, type=type, ng=ng, cp.matrix=cp.matrix)
  Freq=mp[[2]]
  QTL=mp[[3]]
  marker=mp[[4]]
  n.para=ncol(D.matrix)
  Qn=nrow(D.matrix)
  
  L0=sum(log(dnorm(Y,mean=mean(Y),sd=var(Y)^0.5)))
  if(is.null(X)){
    X=matrix(1,N,1)
  }else if(is.vector(X)){
    X=matrix(X,N,1)
  }
  if(is.null(beta)){
    mut=matrix(rep(mean(Y),ncol(X)),ncol(X),1)
  }else if(is.numeric(beta)){
    mut=matrix(rep(beta,ncol(X)),ncol(X),1)
  }else{mut=beta}
  if(is.null(E.vector)){
    Et=rep(0,n.para)
  }else{Et=E.vector}
  if(is.null(variance)){
    sigt=var(Y)
  }else{sigt=variance}
  MUt=t(D.matrix%*%Et%*%matrix(1,1,N))+X%*%mut%*%matrix(1,1,Qn)
  PIup=matrix(0,N,Qn)
  for(j in 1:N){
    P0=c()
    for(i.p in 1:Qn){
      P0[i.p]=Freq[j,i.p]*dnorm(Y[j],mean=MUt[j,i.p],sd=sigt^0.5)
    }
    if(sum(P0)!=0){P0=P0/sum(P0)
    } else {P0=rep(1/g,g)}
    PIup[j,]=P0
  }
  PI=PIup
  
  R=matrix(0,n.para,1)
  M=V=matrix(0,n.para,n.para)
  one=rep(1,N)
  for(m1 in 1:n.para){
    R[m1,]=(t(Y-X%*%mut)%*%PI%*%D.matrix[,m1])/(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))
    for(m2 in 1:n.para){
      V[m1,m2]=one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2])
      if(m1!=m2)
      {M[m1,m2]=(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2]))/
        (one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))}
    }
  }
  Et1=R-M%*%Et
  mut1=solve(t(X)%*%X)%*%t(X)%*%(Y-PI%*%D.matrix%*%Et1)
  sigt1=(t(Y-X%*%mut1)%*%(Y-X%*%mut1)-2*(t(Y-X%*%mut1)%*%PI%*%D.matrix%*%Et1)+(t(Et1)%*%V%*%Et1))/N
  MUt1=t(D.matrix%*%Et1%*%matrix(1,1,N))+X%*%mut1%*%matrix(1,1,Qn)
  Lih10=Lih1=matrix(0,N,Qn)
  for(i.L in 1:Qn){
    Lih10[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt[,i.L],sd=sigt^0.5)
    Lih1[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt1[,i.L],sd=sigt1^0.5)
  }
  L10=sum(log(rowSums(Lih10)))
  L1=sum(log(rowSums(Lih1)))
  time=1
  while(max(abs(c(Et1-Et,sigt1-sigt)))>=conv){
    repeat{
      mut=mut1
      Et=Et1
      sigt=sigt1
      MUt=MUt1
      L10=L1
      PIup=matrix(0,N,Qn)
      for(j in 1:N){
        P0=c()
        for(i.p in 1:Qn){
          P0[i.p]=Freq[j,i.p]*dnorm(Y[j],mean=MUt[j,i.p],sd=sigt^0.5)
        }
        if(sum(P0)!=0){P0=P0/sum(P0)
        } else {P0=rep(1/g,g)}
        PIup[j,]=P0
      }
      PI=PIup
      for(m1 in 1:n.para){
        R[m1,]=(t(Y-X%*%mut)%*%PI%*%D.matrix[,m1])/(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))
        for(m2 in 1:n.para){
          V[m1,m2]=one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2])
          if(m1!=m2)
          {M[m1,m2]=(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2]))/
            (one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))}
        }
      }
      Et1=R-M%*%Et
      mut1=solve(t(X)%*%X)%*%t(X)%*%(Y-PI%*%D.matrix%*%Et1)
      sigt1=(t(Y-X%*%mut1)%*%(Y-X%*%mut1)-2*(t(Y-X%*%mut1)%*%PI%*%D.matrix%*%Et1)+(t(Et1)%*%V%*%Et1))/N
      MUt1=t(D.matrix%*%Et1%*%matrix(1,1,N))+X%*%mut1%*%matrix(1,1,Qn)
      Lih10=Lih1=matrix(0,N,Qn)
      for(i.L in 1:Qn){
        Lih10[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt[,i.L],sd=sigt^0.5)
        Lih1[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt1[,i.L],sd=sigt1^0.5)
      }
      L10=sum(log(rowSums(Lih10)))
      L1=sum(log(rowSums(Lih1)))
      break(max(abs(c(Et1-Et,sigt1-sigt)))<conv)
    }
    if(console==T){
      Ep=round(Et1,5)
      sp=round(sigt1,5)
      print(c("Time",time,"effect",Ep,"var",sp))
    }
    time=time+1
  }
  time=time-1
  LRT=-2*(L0-L1)
  parameter=c(mut1,Et1,sigt1)
  
  Et1=parameter[2:(1+ncol(D.matrix))]
  mut1=parameter[1]
  sigt1=parameter[length(parameter)]
  colnames(PI)=colnames(mp)
  output=list(QTL=QTL,E.vector=Et1,beta=mut1,variance=sigt1,PI.matrix=PI,LRT=LRT,iteration.time=time)
  return(output)
}



QTL.t = function(QTL, marker, geno, D.matrix, ys, yu, tL, tR, cM=T, type="RI", ng=2, 
                 E.vector=NULL, beta=NULL,  X=NULL, variance=NULL, cp.matrix=NULL, conv = 10^-5, console=F){
  maxit = 1000
  nS = length(ys)
  nType = nrow(D.matrix)
  nE = ncol(D.matrix)
  mp=mixprop(QTL, nu, marker, geno, model=2, cM=cM, type=type, ng=ng, cp.matrix=cp.matrix)
  Freq=mp[[2]]
  QTL=mp[[3]]
  marker=mp[[4]]
  Y=c(ys,yu)
  init = c(mean(Y), rep(0, nE), sd(Y))
  negLL0 = function(theta) {
    mu = theta[1]
    sigma = theta[2]
    tauL = (tL - mu)/sigma
    tauR = (tR - mu)/sigma
    U = 1 + pnorm(tauL) - pnorm(tauR)
    LL = sum(log(dnorm(ys, mu, sigma)/U))
    return(-LL)
  }
  temp = optim(c(mean(Y), sd(Y)), negLL0, method = "Nelder-Mead", control=list(maxit = 3000))
  LL0 = -1 * (temp$value)
  
  err = 1
  if(is.null(X)){
    X=matrix(1,nS,1)
    if(length(beta>1)){
      beta=mean(beta)
    }
  }else if(is.vector(X)){
    X=matrix(1,nS,1)
  }
  if(is.null(beta)){
    mu0=matrix(rep(mean(Y),ncol(X)),ncol(X),1)
  }else if(is.numeric(beta)){
    mu0=matrix(rep(beta,ncol(X)),ncol(X),1)
  }else{mut=beta}
  if(is.null(E.vector)){
    E=rep(0,nE)
  }else{E=E.vector}
  if(is.null(variance)){
    sigma=sd(Y)
  }else{sigma=sqrt(variance)}
  mu=t(D.matrix%*%E%*%matrix(1,1,nS))+X%*%mu0%*%matrix(1,1,nType)
  
  tauL = (tL - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
  tauR = (tR - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
  U = 1 + pnorm(tauL) - pnorm(tauR)
  LL = sum(log(rowSums(Freq * dnorm(ys, mu, sigma)/U)))
  Et=matrix(NA,3000,nE)
  Mt=c()
  St=c()
  for (iter in 1:maxit) {
    E0=E
    sigma0=sigma
    mu1=mu0
    Pi = Freq * dnorm(ys, mu, sigma)/U
    Pi[rowSums(Pi)==0,]=1
    Pi = Pi/rowSums(Pi)
    Amy = (dnorm(tauL) - dnorm(tauR))/U
    Bob = (tauL * dnorm(tauL) - tauR * dnorm(tauR))/U
    V = matrix(0, nE, nE)
    for (i in 1:nE) {
      V[i, ] = colSums(Pi) %*% (D.matrix[, i] * D.matrix)
    }
    M = V/diag(V) - diag(1, nE)
    r = (t(ys - X%*%mu0) %*% Pi + sigma * colSums(Pi * Amy)) %*%  D.matrix/colSums(Pi %*% D.matrix^2)
    E = t(r) - M %*% E
    mu0 = solve(t(X)%*%X)%*%t(X)%*%(ys - Pi %*% D.matrix %*% E + sigma * rowSums(Pi * Amy))
    sigma = (t(ys - X%*%mu0) %*% (ys - X%*%mu0) - 2 * t(ys - X%*%mu0) %*% Pi %*%  D.matrix %*% E + t(E) %*% V %*% E)/(nS - sum(Pi * Bob))
    sigma = sqrt(as.vector(sigma))
    mu = t(D.matrix%*%E%*%matrix(1,1,nS))+X%*%mu0%*%matrix(1,1,nType)
    tauL = (tL - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
    tauR = (tR - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
    U = 1 + pnorm(tauL) - pnorm(tauR)
    LL1 = sum(log(rowSums(Freq * dnorm(ys, mu, sigma)/U)))
    err = abs(LL1 - LL)
    LL = LL1
    if(console==T){
      Ep=round(E,5)
      sp=round(sigma^2,5)
      print(c("Time",iter,"effect",Ep,"var",sp))
    }
    Et[iter,]=E
    Mt[iter]=mu0
    St[iter]=sigma
    if (max(abs(c(E-E0,mu0-mu1,sigma-sigma0))) < conv) {
      break
    }
    if(iter==maxit){
      if(ncol(Et)>1){
        E=apply(Et[(maxit-119):maxit,],2,mean)
      } else {E=mean(Et)}
      mu0=mean(Mt[(maxit-119):maxit])
      sigma=mean(St[(maxit-119):maxit])
    }
  }
  parameter = c(mu0, E, sigma^2)
  time=iter
  
  LRT = -2*(LL0 - LL1)
  Et1=parameter[2:(1+ncol(D.matrix))]
  mut1=parameter[1]
  sigt1=parameter[length(parameter)]
  colnames(Pi)=colnames(mp)
  output=list(QTL=QTL,E.vector=Et1,beta=mut1,variance=sigt1,PI.matrix=Pi,LRT=LRT,iteration.time=time)
  return(output)
}




EM.IM2=function(marker,geno,y,yu=NULL,sele.g="n",tL=NULL,tR=NULL,method="EM",type="RI",ng=2,
                D.matrix=NULL,cM=T,speed=1,conv = 10^-5,d.eff=T,LRT.thre=T,simu=1000,lv=0.05,
                detect=T,minQTLdist=15,chart = "chr",console = F){
  if(cM==F){
    marker[,2]=marker[,2]*100
  }
  
  if(is.null(D.matrix)){
    if(type=="BC"){
      D.matrix=matrix(c(0.5,-0.5),2,1)
      row.names(D.matrix)=c(2,1)
      colnames(D.matrix)="a1"
    } else {
      a1=c(1,0,-1)
      d1=c(-0.5,0.5,-0.5)
      if(d.eff==T){
        D.matrix=cbind(a1,d1)
        row.names(D.matrix)=c(2,1,0)
      } else {
        D.matrix=matrix(a1,3,1)
        row.names(D.matrix)=c(2,1,0)
        colnames(D.matrix)="a1"
      }
    }
  }
  
  if(method=="EM"){
    meth=function(QTL,marker,geno,D.matrix,y,yu,tL,tR,type,ng,sele.g,conv){
      EM=EM.MIM2(QTL,marker,geno,D.matrix,y=y,yu=yu,tL=tL,tR=tR,type=type,ng=ng,sele.g=sele.g,conv=conv)
      eff=as.numeric(EM$E.vector)
      mu0=as.numeric(EM$beta)
      sigma=sqrt(as.numeric(EM$variance))
      LRT=EM$LRT
      model=EM$model
      result=list(eff,mu0,sigma,LRT,model)
      return(result)
    }
    Y=y[!is.na(y)]
    Yu=yu[!is.na(yu)]
    if(sele.g=="p"){
      Ya=c(y,yu)
    }else if(sele.g=="t"){
      Ya=y
    }else{
      Ya=y
    }
  } else if (method=="REG"){
    Y=y[!is.na(y)]
    Yu=yu[!is.na(yu)]
    if(sele.g=="p"){
      Ya=c(y,yu)
    }else if(sele.g=="t"){
      Ya=y
    }else{
      Ya=y
    }
    meth=function(QTL,marker,geno,D.matrix,y,yu,tL,tR,type,ng,sele.g,conv){
      if(sele.g=="p"){
        mp=mixprop(QTL, length(Yu), marker, geno, model=1, cM=cM, type=type, ng=ng)[[2]]
      }else if(sele.g=="t"){
        mp=mixprop(QTL, length(Yu), marker, geno, model=2, cM=cM, type=type, ng=ng)[[2]]
      }else{
        mp=Q.make(QTL,marker,geno,type=type,ng=ng)$cp.matrix
      }
      X=mp%*%D.matrix
      fit=lm(Ya~X)
      eff=as.numeric(fit$coefficients[-1])
      mu0=as.numeric(fit$coefficients[1])
      ms=anova(fit)$`Mean Sq`
      sigma=ms[2]^0.5
      
      L0=c()
      L1=c()
      for(k in 1:nrow(mp)){
        L00=c()
        L01=c()
        for(m in 1:nrow(D.matrix)){
          L00[m]=mp[k,m]*dnorm((Ya[k]-mu0)/sigma)
          L01[m]=mp[k,m]*dnorm((Ya[k]-(mu0+D.matrix[m,]%*%eff))/sigma)
        }
        L0[k]=sum(L00)
        L1[k]=sum(L01)
      }
      
      LRT=-2*sum(log(L0[!is.na(L0) & !is.na(L1)]/L1[!is.na(L0) & !is.na(L1)]))
      
      result=list(eff,mu0,sigma,LRT,model="regression interval mapping model")
      return(result)
    }
  }
  
  effect=c()
  cr0=unique(marker[,1])
  for(i in cr0){
    cr=marker[marker[,1]==i,]
    for(j in seq(ceiling(floor(min(cr[,2]))+speed),(max(cr[,2])),speed)){
      QTL=matrix(c(i,j),1,2)
      
      result=meth(QTL,marker,geno,D.matrix,y,yu,tL,tR,type,ng,sele.g,conv)
      eff=result[[1]]
      mu0=result[[2]]
      sigma=result[[3]]
      LRT=result[[4]]
      model=result[[5]]
      
      eff0=c(i,j,eff,LRT)
      effect=rbind(effect,eff0)
      if(console==T){print(c(paste("chr",i),paste(j,"cM"),paste("LRT",LRT)))}
    }
  }
  row.names(effect)=1:nrow(effect)
  colnames(effect)=c("chr","cM",colnames(D.matrix),"LRT")
  effect=data.frame(effect)
  effect[effect[,ncol(effect)]==Inf,3:ncol(effect)]=0
  
  
  LRT.threshold=NULL
  if(LRT.thre==T){
    LRT.threshold=LRTthre(marker,type,ng,length(Ya),25,cM,d.eff,simu,speed,lv,console)
  } else if (class(LRT.thre)=="numeric"){
    LRT.threshold=LRT.thre
  }
  
  
  detect.QTL=NULL
  if(detect==T & class(LRT.threshold)=="numeric"){
    detect.QTL=detectQTL(effect,LRT.threshold,minQTLdist)
  }
  
  if(chart=="chr"){
    for(k in cr0){
      par(mar=c(5,5,4,2))
      par(fig=c(0,1,0.35,1))
      plot(effect[effect$chr==k,]$cM,effect[effect$chr==k,]$LRT,main="",ylab="LRT statistic",
           xlab="",type="l",ylim=c(0,max(c(effect$LRT,LRT.threshold),na.rm = T)),bty="l")
      if(!is.null(LRT.threshold)){abline(h=LRT.threshold,col="red",lwd=2)}
      par(mar=c(2,5,4,2))
      par(fig=c(0,1,0,0.45),new=T)
      plot(effect[effect$chr==k,]$cM,effect[effect$chr==k,]$a1,ylab="effect",main=paste("Chromosome",k),
           xlab="",type="l",ylim=c(min(c(effect$a1,0))*1.4,max(c(effect$a1,0))*1.2),xaxt="n",bty="n")
      abline(h=0,col="black",lwd=2)
    }
  } else if (chart=="all") {
    t=table(effect$chr)
    ncr=cumsum(t)
    par(mar=c(5,5,4,2))
    par(fig=c(0,1,0.3,1))
    plot(effect$LRT,main="",ylab="LRT statistic",xlab="",type="n",
         ylim=c(0,max(c(effect$LRT,LRT.threshold),na.rm = T)),xaxt="n",bty="l")
    k0=0
    for(k in cr0){
      cMk=max(effect[effect$chr==k,]$cM)
      points(effect[effect$chr==k,]$cM+k0,effect[effect$chr==k,]$LRT,main="",ylab="LRT statistic",xlab="",
             type="l",ylim=c(0,max(c(effect$LRT,LRT.threshold),na.rm = T)),xaxt="n",bty="l")
      k0=k0+cMk
      abline(v=k0,col="blue",lwd=2)
    }
    abline(v=k0,col="white",lwd=2)
    axis(1,paste("chr",(1:length(ncr))),at=ncr-t/2)
    if(!is.null(LRT.threshold)){abline(h=LRT.threshold,col="red",lwd=2)}
    
    par(mar=c(2,5,4,2))
    par(fig=c(0,1,0,0.5),new=T)
    plot(effect$a1,main="",ylab="effect",xlab="",type="n",
         ylim=c(min(c(effect$a1,0))*1.2,max(c(effect$a1,0))*1.2),xaxt="n",bty="n")
    k0=0
    for(k in cr0){
      cMk=max(effect[effect$chr==k,]$cM)
      points(effect[effect$chr==k,]$cM+k0,effect[effect$chr==k,]$a1,main="",ylab="effect",xlab="",type="l",
             ylim=c(min(c(effect$a1,0))*1.2,max(c(effect$a1,0))*1.2),xaxt="n",bty="n")
      k0=k0+cMk
      abline(v=k0,col="blue",lwd=2)
    }
    abline(v=k0,col="white",lwd=2)
    abline(h=0,col="black",lwd=2)
  }
  par(mar = c(5.1,4.1,4.1,2.1))
  par(fig=c(0,1,0,1))
  result=list(effect=effect,LRT.threshold=LRT.threshold,detect.QTL=detect.QTL,model=model)
}



#######################################################################################################################

QTL.estim.p=function(QTL, marker, geno, D.matrix, ys, yu, cM=T, type="RI", ng=2, 
                     E.vector=NULL, X=NULL, beta=NULL, variance=NULL, cp.matrix=NULL, EM=T, conv = 10^-5, console=F){
  Y=c(ys,yu)
  N=length(Y)
  nu=length(yu)
  nQTL=nrow(QTL)
  mp=mixprop(QTL, nu, marker, geno, model=1, cM=cM, type=type, ng=ng, cp.matrix=cp.matrix)
  Freq=mp[[2]]
  QTL=mp[[3]]
  marker=mp[[4]]
  n.para=ncol(D.matrix)
  Qn=nrow(D.matrix)
  if(EM==T){
    L0=sum(log(dnorm(Y,mean=mean(Y),sd=var(Y)^0.5)))
    if(is.null(X)){
      X=matrix(1,N,1)
    }else if(is.vector(X)){
      X=matrix(1,N,1)
    }
    if(is.null(beta)){
      mut=matrix(rep(mean(Y),ncol(X)),ncol(X),1)
    }else if(is.numeric(beta)){
      mut=matrix(rep(beta,ncol(X)),ncol(X),1)
    }else{mut=beta}
    if(is.null(E.vector)){
      Et=rep(0,n.para)
    }else{Et=E.vector}
    if(is.null(variance)){
      sigt=var(Y)
    }else{sigt=variance}
    MUt=t(D.matrix%*%Et%*%matrix(1,1,N))+X%*%mut%*%matrix(1,1,Qn)
    PIup=matrix(0,N,Qn)
    for(j in 1:N){
      P0=c()
      for(i.p in 1:Qn){
        P0[i.p]=Freq[j,i.p]*dnorm(Y[j],mean=MUt[j,i.p],sd=sigt^0.5)
      }
      if(sum(P0)!=0){P0=P0/sum(P0)
      } else {P0=rep(1/g,g)}
      PIup[j,]=P0
    }
    PI=PIup
    
    R=matrix(0,n.para,1)
    M=V=matrix(0,n.para,n.para)
    one=rep(1,N)
    for(m1 in 1:n.para){
      R[m1,]=(t(Y-X%*%mut)%*%PI%*%D.matrix[,m1])/(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))
      for(m2 in 1:n.para){
        V[m1,m2]=one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2])
        if(m1!=m2)
        {M[m1,m2]=(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2]))/
          (one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))}
      }
    }
    Et1=R-M%*%Et
    mut1=solve(t(X)%*%X)%*%t(X)%*%(Y-PI%*%D.matrix%*%Et1)
    sigt1=(t(Y-X%*%mut1)%*%(Y-X%*%mut1)-2*(t(Y-X%*%mut1)%*%PI%*%D.matrix%*%Et1)+(t(Et1)%*%V%*%Et1))/N
    MUt1=t(D.matrix%*%Et1%*%matrix(1,1,N))+X%*%mut1%*%matrix(1,1,Qn)
    Lih10=Lih1=matrix(0,N,Qn)
    for(i.L in 1:Qn){
      Lih10[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt[i.L],sd=sigt^0.5)
      Lih1[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt1[i.L],sd=sigt1^0.5)
    }
    L10=sum(log(rowSums(Lih10)))
    L1=sum(log(rowSums(Lih1)))
    time=1
    while(max(abs(c(Et1-Et,sigt1-sigt)))>=conv){
      repeat{
        mut=mut1
        Et=Et1
        sigt=sigt1
        MUt=MUt1
        L10=L1
        PIup=matrix(0,N,Qn)
        for(j in 1:N){
          P0=c()
          for(i.p in 1:Qn){
            P0[i.p]=Freq[j,i.p]*dnorm(Y[j],mean=MUt[j,i.p],sd=sigt^0.5)
          }
          if(sum(P0)!=0){P0=P0/sum(P0)
          } else {P0=rep(1/g,g)}
          PIup[j,]=P0
        }
        PI=PIup
        for(m1 in 1:n.para){
          R[m1,]=(t(Y-X%*%mut)%*%PI%*%D.matrix[,m1])/(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))
          for(m2 in 1:n.para){
            V[m1,m2]=one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2])
            if(m1!=m2)
            {M[m1,m2]=(one%*%PI%*%(D.matrix[,m1]*D.matrix[,m2]))/
              (one%*%PI%*%(D.matrix[,m1]*D.matrix[,m1]))}
          }
        }
        Et1=R-M%*%Et
        mut1=solve(t(X)%*%X)%*%t(X)%*%(Y-PI%*%D.matrix%*%Et1)
        sigt1=(t(Y-X%*%mut1)%*%(Y-X%*%mut1)-2*(t(Y-X%*%mut1)%*%PI%*%D.matrix%*%Et1)+(t(Et1)%*%V%*%Et1))/N
        MUt1=t(D.matrix%*%Et1%*%matrix(1,1,N))+X%*%mut1%*%matrix(1,1,Qn)
        Lih10=Lih1=matrix(0,N,Qn)
        for(i.L in 1:Qn){
          Lih10[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt[i.L],sd=sigt^0.5)
          Lih1[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt1[i.L],sd=sigt1^0.5)
        }
        L10=sum(log(rowSums(Lih10)))
        L1=sum(log(rowSums(Lih1)))
        break(abs(L10-L1)<conv)
      }
      if(console==T){
        Ep=round(Et1,5)
        sp=round(sigt1,5)
        print(c("Time",time,"effect",Ep,"var",sp))
      }
      time=time+1
    }
    m0="EM"
    time=time-1
    LRT=-2*(L0-L1)
    parameter=c(mut1,Et1,sigt1)
  }
  if(EM==F){
    lik0=function(theta){
      mu0=theta[1]
      sig0=theta[2]
      likf0=-sum(log(dnorm(Y,mean=mu0,sd=sig0^0.5)))
      return(likf0)
    }
    L0=optim(c(mean(Y),var(Y)),lik0)
    h0lrt=L0$value
    lik1=function(theta){
      mut=theta[1]
      sigt=theta[n.para+2]
      Et=theta[-(c(1,n.para+2))]
      MUt=c(rep(mut,Qn))+D.matrix%*%matrix(Et,n.para,1)
      Lihood=matrix(0,N,Qn)
      for(i.L in 1:Qn){
        Lihood[,i.L]=Freq[,i.L]*dnorm(Y,mean=MUt[i.L],sd=sigt^0.5)
      }
      likf1=-sum(log(rowSums(Lihood)))
      return(likf1)
    }
    L1=optim(c(mean(Y),rep(0,n.para),var(Y)),lik1,control=list(maxit=3000))
    parameter=L1$par
    Et1=parameter[1:ncol(D.matrix)]
    mut1=parameter[ncol(D.matrix)+1]
    sigt1=parameter[ncol(D.matrix)+2]
    m0="optim"
    time=NULL
    PI=NULL
    h1lrt=L1$value
    LRT=2*(h0lrt-h1lrt)
  }
  Et1=parameter[2:(1+ncol(D.matrix))]
  mut1=parameter[1]
  sigt1=parameter[length(parameter)]
  output=list(E.vector=Et1,beta=mut1,varence=sigt1,method=m0,PI.matrix=PI,time=time)
  return(output)
}


QTL.estim.t = function(QTL, marker, geno, D.matrix, ys, tL, tR, cM=T, type="RI", ng=2, 
                       E.vector=NULL, beta=NULL,  X=NULL, variance=NULL, cp.matrix=NULL, EM=T, conv = 10^-5, console=F){
  maxit = 3000
  nS = length(ys)
  nType = nrow(D.matrix)
  nE = ncol(D.matrix)
  mp=mixprop(QTL, nu, marker, geno, model=2, cM=cM, type=type, ng=ng, cp.matrix=cp.matrix)
  Freq=mp[[2]]
  QTL=mp[[3]]
  marker=mp[[4]]
  init = c(mean(ys), rep(0, nE), sd(ys))
  negLL0 = function(theta) {
    mu = theta[1]
    sigma = theta[2]
    tauL = (tL - mu)/sigma
    tauR = (tR - mu)/sigma
    U = 1 + pnorm(tauL) - pnorm(tauR)
    LL = sum(log(dnorm(ys, mu, sigma)/U))
    return(-LL)
  }
  temp = optim(c(mean(ys), sd(ys)), negLL0, method = "Nelder-Mead", control=list(maxit = maxit))
  LL0 = -1 * (temp$value)
  if (EM==T) {
    err = 1
    if(is.null(X)){
      X=matrix(1,nS,1)
      if(length(beta>1)){
        beta=mean(beta)
      }
    }else if(is.vector(X)){
      X=matrix(1,nS,1)
    }
    if(is.null(beta)){
      mu0=matrix(rep(mean(ys),ncol(X)),ncol(X),1)
    }else if(is.numeric(beta)){
      mu0=matrix(rep(beta,ncol(X)),ncol(X),1)
    }else{mut=beta}
    if(is.null(E.vector)){
      E=rep(0,nE)
    }else{E=E.vector}
    if(is.null(variance)){
      sigma=sd(ys)
    }else{sigma=sqrt(variance)}
    mu=t(D.matrix%*%E%*%matrix(1,1,nS))+X%*%mu0%*%matrix(1,1,nType)
    
    tauL = (tL - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
    tauR = (tR - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
    U = 1 + pnorm(tauL) - pnorm(tauR)
    LL = sum(log(rowSums(Freq * dnorm(ys, mu, sigma)/U)))
    for (iter in 1:maxit) {
      E0=E
      sigma0=sigma
      mu1=mu0
      Pi = Freq * dnorm(ys, mu, sigma)/U
      Pi[rowSums(Pi)==0,]=1
      Pi = Pi/rowSums(Pi)
      Amy = (dnorm(tauL) - dnorm(tauR))/U
      Bob = (tauL * dnorm(tauL) - tauR * dnorm(tauR))/U
      V = matrix(0, nE, nE)
      for (i in 1:nE) {
        V[i, ] = colSums(Pi) %*% (D.matrix[, i] * D.matrix)
      }
      M = V/diag(V) - diag(1, nE)
      r = (t(ys - X%*%mu0) %*% Pi + sigma * colSums(Pi * Amy)) %*%  D.matrix/colSums(Pi %*% D.matrix^2)
      E = t(r) - M %*% E
      mu0 = solve(t(X)%*%X)%*%t(X)%*%(ys - Pi %*% D.matrix %*% E + sigma * rowSums(Pi * Amy))
      sigma = (t(ys - X%*%mu0) %*% (ys - X%*%mu0) - 2 * t(ys - X%*%mu0) %*% Pi %*%  D.matrix %*% E + t(E) %*% V %*% E)/(nS - sum(Pi * Bob))
      sigma = sqrt(as.vector(sigma))
      mu = t(D.matrix%*%E%*%matrix(1,1,nS))+X%*%mu0%*%matrix(1,1,nType)
      tauL = (tL - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
      tauR = (tR - X%*%mu0%*%matrix(1,1,nType) - t(D.matrix%*%E%*%matrix(1,1,nS)))/sigma
      U = 1 + pnorm(tauL) - pnorm(tauR)
      LL1 = sum(log(rowSums(Freq * dnorm(ys, mu, sigma)/U)))
      err = abs(LL1 - LL)
      LL = LL1
      if(console==T){
        Ep=round(E,5)
        sp=round(sigma^2,5)
        print(c("Time",iter,"effect",Ep,"var",sp))
      }
      if (max(abs(c(E-E0,mu0-mu1,sigma-sigma0))) < conv) {
        break
      }
    }
    parameter = c(mu0, E, sigma^2)
    m0="EM"
    time=iter
  }
  if (EM==F) {
    negLL1 = function(theta) {
      mu0 = theta[1]
      E = theta[1 + 1:nE]
      sigma = tail(theta, 1)
      mu = mu0 + matrix(D.matrix %*% E, nS, nType, byrow = TRUE)
      tauL = (tL - mu0 - D.matrix %*% E)/sigma
      tauR = (tR - mu0 - D.matrix %*% E)/sigma
      U = matrix(1 + pnorm(tauL) - pnorm(tauR), nS, nType, byrow = TRUE)
      LL = sum(log(rowSums(Freq * (dnorm(ys, mu, sigma)/U))))
      return(-LL)
    }
    temp = optim(init, negLL1, method = "Nelder-Mead", control=list(maxit=maxit, abstol=conv))
    LL1 = -1 * (temp$value)
    parameter = temp$par
    parameter[nE + 2] = parameter[nE + 2]^2
    m0="optim"
    time=NULL
    Pi=NULL
  }
  LRT = -2*(LL0 - LL1)
  Et1=parameter[2:(1+ncol(D.matrix))]
  mut1=parameter[1]
  sigt1=parameter[length(parameter)]
  output=list(E.vector=Et1,beta=mut1,varence=sigt1,method=m0,PI.matrix=Pi,time=time)
  return(output)
}