Nothing
#for now, history is given as matrix
#(1/ncopies, ..., 1); Must be length equal to the number of copies of chromosome for the history
#later will have function that takes encoding of history and creates this vector
#if exactAllele=TRUE, then assume that x/m is exactly the true allele frequency P -- not implemented yet
#... corresponds to arguments passed to the internal fitting functions (e.g. Bayes method of approximation)
###Should code some EM steps to get good init point...
##Version 1.0.2 add ids for mutations; made returnAssignments and returnData the same thing -- one data.frame
eventTiming<-function(x, m, history, totalCopy, method=c("fullMLE","partialMLE","Bayes"),type=c("gain","CNLOH"),
seqError=0, bootstrapCI=NULL, B=if(method=="Bayes") 10000 else 500,CILevel=0.95, normCont=0,verbose=TRUE,
returnAssignments=FALSE,coverageCutoff=1,minMutations=10,init=NULL,maxiter=100, tol=0.0001,mutationId=1:length(x),...){
method<-match.arg(method)
doCI<-!is.null(bootstrapCI)
type<-match.arg(type)
nMuts<-length(x)
if(length(m)!=nMuts) stop("x and m must be of same length")
if(length(mutationId)!=length(unique(mutationId)) & verbose) warning("mutationId not unique values; using index as id instead.")
if(length(mutationId)!=nMuts & verbose) warning("mutationId an invalid length; using index as id instead.")
if(returnAssignments) returnData<-TRUE
else returnData<-FALSE
A<-history #this will be a function later
nCopies<-totalCopy
if(is.null(dim(A))) stop("'history' should be a matrix of size (nEvents +1) x (nEvents +1)")
if(ncol(A)!=nrow(A)) stop("'history' should be a square matrix (hint: do not include the allele frequency 1 except for CNLOH events)")
K<-ncol(A)-1
if(length(normCont)!=1 || normCont>1 || normCont<0) stop("normCont must be given and should be single value between 0 and 1.")
#get final allele frequencies
possAlleles<-allAF(nCopies,normCont=normCont,type="mutation")[[1]]
if(type=="gain") possAlleles<-head(possAlleles,-1) #can't have the allele frequency 1
if(nrow(A)!=length(possAlleles)) stop(length(possAlleles),"possible alleles for this scenario, but history only has",nrow(A),"rows.")
possAlleles<-errorAF(possAlleles,seqError=seqError)
if(coverageCutoff<1) stop("coverageCutoff must be at least 1")
nBelowCutoff<-length(which(m<coverageCutoff))
nZero<-length(which(m==0))
if(any(m<coverageCutoff)){
if(verbose) warning(nBelowCutoff, "mutations present below cutoff",coverageCutoff,", will be ignored")
if(verbose & nZero>0) warning(nZero, "mutations have no reads at all")
x<-x[which(m>=coverageCutoff)]
mutationId<-mutationId[which(m>=coverageCutoff)]
m<-m[which(m>=coverageCutoff)]
}
summaryTable<-cbind(nMuts,length(m),nZero,nBelowCutoff)
colnames(summaryTable)<-c("Total Mutations","Qualify for fitting","Zero Read Coverage","Coverage Below Cutoff")
summaryTable<-t(summaryTable)
colnames(summaryTable)<-"Number"
call<-list(alleleSet=possAlleles,
history=history,
totalCopy=totalCopy,
type=type,
method=method,
seqError=seqError,
normCont=normCont,
coverageCutoff=coverageCutoff,
minMutations=minMutations,
B=B,
init=init,
maxiter=maxiter,
tol=tol)
if(returnData) inputData=data.frame(mutationId=mutationId,x=x,m=m)
failReason<-NULL
success<-TRUE
if(length(m)<minMutations){
failReason<-paste(failReason,"not enough mutations above cutoff that satisfy coverage criteria",sep=",")
if(verbose) warning("less than ",minMutations," mutated locations meet coverage criteria; no estimate of pi will be calculated")
success<-FALSE
}
# if(all(x==m)){
# failReason<-paste(failReason,"all x equal to m",sep=",")
# success<-FALSE
# }
dummyOutput<-list("pi"=rep(NA,length=ncol(A)),alleleSet=possAlleles,optimDetails=list(nIter=NA,initial=NA,pathOfQ=NA),perLocationProb=NA,q=NA)
names(dummyOutput$pi)<-paste("Stage",0:K,sep="")
ci<-matrix(rep(NA,length=2*length(dummyOutput$pi)),ncol=2)
row.names(ci)<-names(dummyOutput$pi)
dummyOutput$piCI<-ci
if(!is.null(failReason)){
output<-dummyOutput
}
else{
rankA<-qr(A)$rank
piId <- rankA==ncol(A)
piEst<-rep(NA,ncol(A))
if(piId){
if(method%in%c("fullMLE","Bayes")){
### Estimate the q vector:
output<-try(.estimateQ(x,m,possAlleles,history=history,init=init,maxiter=maxiter, tol=tol))
if(output$optimDetails$nIter==maxiter){
failReason<-paste(failReason,"hit maximum number of iterations",sep=",")
success<-FALSE
}
}
if(method=="partialMLE" ){
mleAlleles<-mleAF(x,m,totalCopy=nCopies,seqError=seqError,normCont=normCont,maxCopy=if(type=="gain") totalCopy-1 else totalCopy)
#q<-mleAlleles$alleleSet$Freq/sum(mleAlleles$alleleSet$Freq)
output<-try(.estimateQ(x,m,possAlleles,alleleFreq=mleAlleles$alleleSet$Freq,history=history,init=init,maxiter=maxiter, tol=tol))
# output<-list(q=q,perLocationProb=mleAlleles$perLocationProb,alleleSet=possAlleles,optimDetails=list(nIter=NA,initial=NA,pathOfQ=NA))
}
if(!inherits(output, "try-error")){
piEst<-solve(A)%*%output$q
piEst<-as.vector(piEst)
piEst<-piEst/sum(piEst)
if(method=="Bayes"){
initBayes<-piEst
#initBayes<-seq(1,length(possAlleles))/length(possAlleles)
bayesout<-.bayesEstimate(x,m,alleleSet=possAlleles,history=A,init=initBayes,nsim=B,CILevel=CILevel,...)
piEst<-bayesout$pi
q<-A%*%piEst
q<-q/sum(q)
output<-list(piCI=bayesout$piCI,q=q,perLocationProb=NA,alleleSet=possAlleles,optimDetails=list(nIter=NA,initial=rbind(initBayes,bayesout$mode),pathOfQ=NA))
}
}
else{
if(method=="Bayes"){
initBayes<-seq(1,length(possAlleles))/length(possAlleles)
bayesout<-.bayesEstimate(x,m,alleleSet=possAlleles,history=A,init=initBayes,nsim=B,...)
piEst<-bayesout$pi
q<-A%*%piEst
q<-q/sum(q)
output<-list(piCI=bayesout$piCI,q=q,perLocationProb=NA,alleleSet=possAlleles,optimDetails=list(nIter=NA,initial=rbind(initBayes,bayesout$mode),pathOfQ=NA))
}
else{
output<-dummyOutput
failReason<-paste(failReason,"error returned in estimating q",sep=",")
success<-FALSE
}
}
if(method=="Bayes") row.names(output$optimDetails$initial)<-c("mle","postMode")
}
else{
output<-dummyOutput
failReason<-paste(failReason,paste("history matrix not invertible, rank=",rankA,sep=""),sep=",")
success<-FALSE
if(verbose) warning(paste("history matrix not invertible, rank=",rankA,sep=""))
}
names(piEst)<-paste("Stage",0:K,sep="")
output$pi<-piEst
if(!returnAssignments) output<-output[-grep("perLocationProb",names(output))]
else{ output[["perLocationProb"]]<-data.frame(inputData,output[["perLocationProb"]],check.names=FALSE)}
#do bootstrap CI, or format the Bayesian ones; makes sure the elements returned are in right order
if(doCI & method!="Bayes"){
if(!any(is.na(output$pi))){
bootsample<-bootstrapEventTiming(call=call,B=B,x=x,m=m,type=bootstrapCI,pi=output$pi)
ci<-t(apply(bootsample,2,quantile,c((1-CILevel)/2,1-(1-CILevel)/2)))
}
else {
ci<-matrix(rep(NA,length=2*length(output$pi)),ncol=2)
row.names(ci)<-names(output$pi)
}
#colnames(ci)<-c("lower","upper")
output$piCI<-ci
output<-output[c("pi","piCI","q","perLocationProb","optimDetails")]
}
else{
if(method!="Bayes") output<-output[c("pi","q","perLocationProb","optimDetails")]
else{
if(!any(is.na(output$pi))){ row.names(output$piCI)<-names(output$pi)}
output<-output[c("pi","piCI","q","perLocationProb","optimDetails")]
}
}
}
return(c(output,list(summaryTable=summaryTable,success=success,failReason=if(!is.null(failReason)) gsub("^,","",failReason) else failReason,call=call)))
}
#simple EM algorithm to calculate \hat{q}, the proportions of the multinomial, constrained
#also gives probabilities of assignments to different alleles.
#if alleleFreq not null, then means know the P[i], and alleleFreq is the number of P[i] equal to each allele
#xGreaterZero determines whether the likelihood will account for the fact that only observe X[i]>0
#useGradient determines whether to use the gradient I calculated, or set gr=NULL
.estimateQ<-function(x,m,alleleSet,alleleFreq=NULL,history,init=NULL,maxiter=100, tol=0.0001,xGreaterZero=TRUE, useGradient=TRUE){
N<-length(x)
possAlleles<-alleleSet
lengthOfTheta<-length(possAlleles)-1
if(nrow(history)!=length(possAlleles)) stop("length of alleleSet must equal the number of rows of history.")
##check if it is one of the easy ones:
Avec<-as.vector(history)
easyType<-"no"
if(identical(Avec,c(0, 1, 2, 0))) easyType<-"CNLOH"
if(identical(Avec,c(1 ,1, 3, 0))) easyType<-"SingleGain"
Ainv<-solve(history)
if(is.null(init)){
init<-history%*%rep(1/length(possAlleles),length=length(possAlleles))
init<-as.vector(init)/sum(init)
# init<-c(.2,.8) #just for testing
}
##Break it into 2 because useful to have a function that assigns most likely allele to each
EStep<-function(q){ ####Calculate P(Pi=a | Xi,q) each i and a
#P(Pi=a | Xi)=P(Xi|Pi=a)P(Pi=a)/sum(P(Xi|Pi=a)P(Pi=a))
#calculate the log of the probability (numerator)
logPPerAllele<-mapply(possAlleles,q,FUN=function(a,qa){
dbinom(x=x, size=m, prob=a,log=TRUE)+log(qa)
}) #returns a N x numb.Alleles Matrix; needs to be normalized
#-----
#calculate log of sum of probabilities (denominator)
#-----
#note if individual probabilities very small, then denominator sums to zero, so divide by zero -- problem!
#use log sum exponential trick:
# log(e[theta1]+e[theta2]+...)=m+log(e[theta1-m]+e[theta2-m]+...), where m=max(theta[i])
#here theta1=log[P(Xi|Pi=a)]+log[P(Pi=a)]
#-----
rowMax<-apply(logPPerAllele,1,max) #max of theta[i]
logPPerAlleleMinusMax<-sweep(logPPerAllele,1,rowMax,"-") #subtract off m from each row
logdenom<-rowMax+log(apply(exp(logPPerAlleleMinusMax),1,sum))
#-----
#calculate log of sum of probabilities (denominator)
#-----
#now overall prob=exp[num]/exp[denom]=exp[num-denom]
#-----
logPPerAllele<-sweep(logPPerAllele,1,logdenom,"-")
# #-----
# #Old simplistic code:
# #-----
# PPerAllele<-mapply(possAlleles,q,FUN=function(a,qa){
# dbinom(x=x, size=m, prob=a,log=FALSE)*(qa)
# }) #returns a N x numb.Alleles Matrix; needs to be normalized
# PPerAllele<- prop.table(PPerAllele,1) #divide by the sum of each row
#compare the two:
# cbind(exp(logPPerAllele),PPerAllele)
return(exp(logPPerAllele))
}
###Constraints for the q
#The feasible region is defined by ui %*% theta - ci >= 0.
if(easyType=="no"){
uiL1<-diag(rep(-1,lengthOfTheta),nrow=lengthOfTheta,ncol=lengthOfTheta)
ciL1<-rep(1,lengthOfTheta)
uiGr0<-diag(rep(1,lengthOfTheta),nrow=lengthOfTheta,ncol=lengthOfTheta)
ciGr0<-rep(0,lengthOfTheta)
uiA<-Ainv%*%(rbind(diag(lengthOfTheta),rep(-1,lengthOfTheta)))
ciA<- - Ainv%*%c(rep(0,lengthOfTheta),1)
ui<-rbind(uiL1,uiGr0,uiA)
ci<- c(-ciL1,-ciGr0,ciA)
}
#make matrix of (1-a[j])^m[i]
MStep<-function(PMat=NULL,qinit){ #estimate of q constrained so that solve(history)%*%q is all positive
if(is.null(alleleFreq)) Na<-colSums(PMat) #from E-Step, the estimated values from the multinomial...basically just the maximum likelihood estimation if knew Pi exactly
else Na<-alleleFreq
if(easyType %in% c("CNLOH","SingleGain")){
#print("easytype")
if(easyType=="CNLOH"){
a<-0
b<-1
}
if(easyType=="SingleGain"){
a<-1/2
b<-1
}
qout<-Na[1]/sum(Na)
if(qout>b) qout<-b
if(qout<a) qout<-a
qout<-c(qout,1-qout)
}
else{
f<-function(q){
qS<-1-sum(q)
#-dmultinom(Na, prob=c(q,qS), log = TRUE) #just to check; should be equivalent, but faster to not calculated normalizing constant
ll<- sum(Na*log(c(q,qS)))
if(xGreaterZero){
#missZero[i]= 1- sum_j{(1-a[j])^m[i] q[j]}
missZero<-1-rowSums(mapply(possAlleles,c(q,qS),FUN=function(a,qa){
(1-a)^m*qa
})#returns a N x numb.Alleles Matrix; needs to be normalized
) #sum over all j, so missZero is vector of length i single value for each i
ll<-ll-sum(missZero)
}
return( - ll) #because minimizing
}
#used in the gradient
#calculate (1-a[S])^{m[i]} - (1-a[j])^{m[i]}
#returns a N x S-1 Matrix; needs to be normalized
aS<-tail(possAlleles,1)
# #bug reported
# num<-mapply(head(possAlleles,-1),q,FUN=function(a,qa){
# (1-aS)^m - (1-a)^m
# })
#fix with
num<-sapply(head(possAlleles,-1),FUN=function(a,qa){
(1-aS)^m - (1-a)^m
})
gr<-function(q){
qS<-1-sum(q)
NS<-tail(Na,1)
g<-(head(Na,-1)/q - NS/qS) #was - (head(Na,-1)/q - NS/qS*q) #which is error; why got wrong answer.
if(xGreaterZero){
#missZero[i]= 1- sum_j{(1-a[j])^m[i] q[j]}
missZero<-1-rowSums(mapply(possAlleles,c(q,qS),FUN=function(a,qa){
(1-a)^m*qa
})#returns a N x numb.Alleles Matrix; needs to be normalized
) #sum over all j, so missZero is vector of length i single value for each i
g<-g-colSums(sweep(num,1,missZero,"/"))
}
return(-g) #because minimizing
}
theta<-head(qinit,-1)
if(lengthOfTheta==1){
upper<-min((ci/ui)[ui<0])
lower<-max((ci/ui)[ui>0])
out<-optimize(f=f,interval=c(lower,upper),tol=.Machine$double.eps)
qout<-c(out$minimum,1-out$minimum)
}
else{
out<-constrOptim(theta=theta,f = f, grad=if(useGradient) gr else NULL, ui=ui, ci=ci) #when use gr, get wrong answer -- returns me to init, doesn't iterate, etc. WHY????
qout<-c(out$par,1-sum(out$par))
}
}
names(qout)<-names(possAlleles)
return(qout)
}
#print("gr uses - (head(Na,-1)/q - NS/qS) ")
TotalStep<-function(q){
MStep(EStep(q),qinit=q)
}
if(is.null(alleleFreq)){
qOld <- init;
qNew<-TotalStep(qOld)
allQ<-cbind(qOld,qNew)
nIter<-1
#convMessages<-c(firstOut$convergence)
while(sum(abs(qNew - qOld)) >= tol & nIter <= maxiter){
qOld<-qNew
pMat<-EStep(qOld)
qNew<-MStep(pMat,qinit=qOld)
#qNew<-c(mstepOut$par,1-sum(mstepOut$par))
#convMessages<-c(convMessages,mstepOut$convergence)
nIter <- nIter + 1
allQ<-cbind(allQ,qNew)
#if(nIter %% 20 ==0) print(sprintf("finished iteration %s", nIter))
}
pMat<-EStep(qNew)
}
else{
qNew<-MStep(qinit=init)
pMat<-NA
nIter<-NA
allQ<-NA
}
names(qNew)<-names(possAlleles)
return(list(q=qNew,perLocationProb=pMat,alleleSet=possAlleles,optimDetails=list(nIter=nIter,initial=init,pathOfQ=allQ)))
}
# > x$pi
# Stage0 Stage1
# 0.1096573 0.8903427
# > x$q
# 1/2 2/2
# 0.96206078 0.03793922
#the Learn Bayes does weird things with the error
#I copy them here to try to fix it; If I comment out the warn options, then will get wanring about optim not appropriate for 1-dim.
.laplace<-function (logpost, mode, ...)
{
currWarn<-options()$warn
#print(currWarn)
options(warn = -1)
fit = optim(mode, logpost, gr = NULL, ..., hessian = TRUE,
control = list(fnscale = -1))
#options(warn = 0)
options(warn = currWarn)
mode = fit$par
h = -solve(fit$hessian)
p = length(mode)
int = p/2 * log(2 * pi) + 0.5 * log(det(h)) + logpost(mode, ...)
stuff = list(mode = mode, var = h, int = int, converge = fit$convergence == 0)
return(stuff)
}
.sir<-function (logf, tpar, n)
{
k = length(tpar$m)
theta = LearnBayes::rmt(n, mean = c(tpar$m), S = tpar$var, df = tpar$df)
lf = matrix(0, c(dim(theta)[1], 1)) #really just n x 1 bc dim(theta)[1]=n
for (j in 1:dim(theta)[1]) lf[j] = logf(theta[j, ])
lp = LearnBayes::dmt(theta, mean = c(tpar$m), S = tpar$var, df = tpar$df,
log = TRUE)
md = max(lf - lp)
wt = exp(lf - lp - md)
probs = wt/sum(wt)
indices = sample(1:n, size = n, prob = probs, replace = TRUE)
if (k > 1)
theta = theta[indices, ]
else theta = theta[indices]
return(theta)
}
.rejectsampling<-function (logf, tpar, dmax, n)
{
d = length(tpar$m)
theta = LearnBayes::rmt(n, mean = c(tpar$m), S = tpar$var, df = tpar$df)
lf = matrix(0, c(dim(theta)[1], 1))
for (j in 1:dim(theta)[1]) lf[j] = logf(theta[j, ])
lg = LearnBayes::dmt(theta, mean = c(tpar$m), S = tpar$var, df = tpar$df,
log = TRUE)
if (d == 1) {
prob = exp(c(lf) - lg - dmax)
return(theta[runif(n) < prob])
}
else {
prob = exp(lf - lg - dmax)
return(theta[runif(n) < prob, ])
}
}
# data(mutData)
# ACNLOH<-matrix(c(1,3,1,0),ncol=2,nrow=2,byrow=TRUE)
# onlyMuts<-subset(mutData,is.na(rsID) & position <= 1.8E7)
# onlyMuts$t_depth<-onlyMuts$t_ref_count+onlyMuts$t_alt_count
# x<-eventTiming(x=onlyMuts$t_alt_count,m=onlyMuts$t_depth,history=ACNLOH,totalCopy=2,type="CNLOH",normCont=0.22)
#outnew<-apply(combSeq,1,function(th){.logPosterior(th,x=xout$nMutAllele,m=xout$nReads,alleleSet= eventOut$call$alleleSet,history=AMat,jacob="new")})
# > length(outnew)
# [1] 10000
# > setwd('/Users/epurdom/Documents/Sequencing/CancerTiming/Simulations')
# > outold<-sapply(xseq,function(th){.logPosterior(th,x=xout$nMutAllele,m=xout$nReads,alleleSet= eventOut$call$alleleSet,history=AMat,jacob="old")})
# xseq<-seq(-10,10,length=100)
# outnew<-sapply(xseq,function(th){.logPosterior(th,x=onlyMuts$t_alt_count,m=onlyMuts$t_depth,alleleSet= x$call$alleleSet,history=ACNLOH,jacob="new")})
# outold<-sapply(xseq,function(th){.logPosterior(th,x=onlyMuts$t_alt_count,m=onlyMuts$t_depth,alleleSet= x$call$alleleSet,history=ACNLOH,jacob="old")})
#run into problems when get close to the boundary --
#conversion between theta and pi is tenuous
#replace log(c+sum(a*exp(theta))) with the following
.mylogSumExp<-function(x,c=1,a=1){
if(length(a)==1) a<-rep(a,length(x))
# if(c==0 || any(exp(x)>1e10)){
# #ignore the +c term in this case and just calculate log( sum_i a_i exp(x_i))
# #then log( sum_i a_iexp(x_i))=m + log(sum_i a_i exp(x_i-m)), where m = max(x_i) http://lingpipe-blog.com/2009/06/25/log-sum-of-exponentials/
# }
#else log(c+sum(a*exp(x))) #this can always be written as c=0 if add term to the x...
if(c!=0){
x<-c(x,log(c))
a<-c(a,1)
}
m<-max(x)
return(m+log(sum(a*exp(x-m))))
}
#always return on the theta scale, but can calculate either log density of theta or pi, depending on choice of value of 'parameter'
.logPosterior<-function(theta,x, m, alleleSet, history, alpha = 1,parameter=c("theta","pi")){
parameter<-match.arg(parameter)
K<-length(theta)
if(length(alleleSet)!=K+1) stop("invalid values for theta or alleleSet")
##############
##Parts of calculating the posterior
pOfPi <- function(theta) {
# (alpha - 1) * (sum(theta - log(1 + sum(exp(theta))))) #change 9/26, but should be equivalent...
(alpha - 1) * (sum(theta) - (K+1)*.mylogSumExp(theta) )
}
pOfX<-function(theta){ #changed 9/30/2013 so that can handle large theta:
PPerX<-mapply(x,m,FUN=function(dat,size){
p<-dbinom(x=dat, size=size, prob=alleleSet) #vector of probabilities (per allele) of seeing data [i] with that allele
#p'A
Avec<-apply(history,2,function(x){sum(x*p)}) #avec[j]=sum_a P_a(X) A[a,j] for each stage j
.mylogSumExp(c(theta,log(1)),c=0,a=Avec) # log( sum_j avec[j]exp(theta[j]) + avec[K] )
})
if(length(PPerX)!=length(x)) stop("coding error")
#Nlog(sum(S[j]*v[j]))
S<-2+0:(K) # this is sequential amplification; could adjust here to be for more complicated
denomTerm<-length(x)*.mylogSumExp(c(theta,log(1)),c=0,a=S) #log( sum_j (2+j)*exp(theta_j) + (2+K))
return(sum(PPerX)-denomTerm)
}
chgOfVar<-function(theta){ #new version 9/26/2013; doesn't affect unless dim(theta)>1
K<-length(theta)
#sum(theta)-(K+1)*log(1+sum(exp(theta)))
sum(theta)-(K+1)*.mylogSumExp(theta)
}
if(parameter=="theta") return(pOfPi(theta) + chgOfVar(theta) + pOfX(theta))
else return(pOfPi(theta) + pOfX(theta))
}
.bayesEstimate<-function (x, m, alleleSet, history, alpha = 1, bayesApproxMethod = c("sir", "inv", "rej","exact"), init, tdf = 4, nsim = 10000, CILevel = 0.95)
{
#change of variables so that now theta[j]=log(pi[j-1]/pi[K])
#assume alpha a single scalar
K<-length(init)-1
method <- match.arg(bayesApproxMethod)
if(method=="inv" & K>1) stop("direct inversion method is only implemented for K=1")
if (!method %in% c("sir","inv"))
stop("Other methods than 'sir' are not yet operational")
# require(LearnBayes) #copied the functions into the code so I could have better control of them, but still need some helper functions...
pi2theta <- function(piV) {
head(log(piV/tail(piV, 1)), -1)
}
theta2pi <- function(theta, returnFullPi=TRUE) {
logPiV<-theta-.mylogSumExp(theta,c=1)
piV<-exp(logPiV)
if(returnFullPi) return(c(piV, 1 - sum(piV)))
else return(piV)
}
if (method == "inv") { #if 1-dimensional, can approximate it easily with a grid of the possibilities
logPiPosterior<-function(theta){.logPosterior(theta,x=x,m=m,alleleSet=alleleSet,history=history,alpha=alpha,parameter="pi")}
piVec<-seq(1e-5,1-1e-5,length=1000)
piVec<-cbind(piVec,1-piVec)
thetaVec<-apply(piVec,1,pi2theta)
#store the values on the log scale so as to make it numerically stable.
logpdensity<- sapply(thetaVec,logPiPosterior)
logNormFac<-.mylogSumExp(logpdensity,c=0,a=1)
logpdensity<-logpdensity-logNormFac
logcumSum<-sapply(1:length(logpdensity),function(ii){
.mylogSumExp(logpdensity[1:ii],c=0)
})
piV<-c(piVec[,1],1)
cumSum<-c(0,exp(logcumSum),1)
#cdf<-stepfun(piV,cumSum,right=TRUE) #f(x[i])=y[i]; otherwise f(x[i])=y[i-1]; #don't need this, actually
meanPi0<-exp(.mylogSumExp(logpdensity,a=piVec[,1],c=0))
upperCred<-uniroot(stepfun(piV,cumSum-(1 - (1 - CILevel)/2),right=TRUE),interval=c(0,1))$root
lowerCred<-uniroot(stepfun(piV,cumSum-(1-CILevel)/2,right=TRUE),interval=c(0,1))$root
ciMat<-rbind(c(lowerCred,upperCred),c(1-upperCred,1-lowerCred))
colnames(ciMat)<-c("2.5%" , "97.5%")
return(list(pi = c(meanPi0,1-meanPi0), piCI = ciMat, modeFit = c(NA,NA),
nsim = nsim)) #return the function that creates the logposterior for checking purposes.
}
else{
logPosterior<-function(theta){.logPosterior(theta,x=x,m=m,alleleSet=alleleSet,history=history,alpha=alpha)}
fit <- try(.laplace(logPosterior, pi2theta(init)))
if (inherits(fit, "try-error")) {
fit <- try(.laplace(logPosterior, pi2theta(rep(1, length(init))/length(init))))
}
tpar <- list(m = fit$mode, var = 2 * fit$var, df = tdf)
if (method == "sir") {
theta.s <- .sir(logPosterior, tpar, nsim)
}
if (method == "rej") {
Tdiff <- function(theta, datapar) {
data <- datapar$data
fit <- datapar$fit
d <- logPosterior(theta, x) - LearnBayes::dmt(theta, mean = fit$mode,
S = 2 * fit$var, df = tdf, log = TRUE)
return(d)
}
maxD <- .laplace(Tdiff, fit$mode, list(data = x, fit = fit))
d <- Tdiff(maxD$mode, list(data = x, fit = fit))
theta.s <- .rejectsampling(logPosterior, tpar, d, nsim)
}
if (is.null(dim(theta.s))) {
pis <- sapply(theta.s, theta2pi)
}
else {
pis <- do.call(cbind, lapply(1:nrow(theta.s), function(i) {
theta2pi(theta.s[i, ])
}))
}
out <- t(apply(pis, 1, quantile, prob = c((1 - CILevel)/2,
0.5, 1 - (1 - CILevel)/2)))
meanOut<- apply(pis, 1, mean)
return(list(pi = meanOut, piCI = out[, c(1, 3)], modeFit = theta2pi(fit$mode),
nsim = nsim)) #return the function that creates the logposterior for checking purposes.
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.