17/11/22 10:27:31.38 WHgLvp/3.net
ab2mv<-function(a,b){
m<-a/(a+b)
v<-m*(1-m)/(a+b+1)
mv<-c(m,v)
return(mv)
}
mv2ab<-function(m,v){
a=(-m^3+m^2-m*v)/v
b=(m^3-2*m^2+m*v+m-v)/v
ab<-c(a,b)
return(ab)
}
HDCI <- function(PMF,cl=0.95){ # Highest Density Confidence Interval
PDF=PMF/sum(PMF)
rsPDF=rev(sort(PDF))
min.density=rsPDF[min(which(cumsum(rsPDF)>=cl))]
index=which(PDF>=min.density)
data.frame(lower.idx=round(min(index)),upper.idx=round(max(index)),actual.CI=sum(PDF[index]))
}