update_pseudo_n_PR: sample alpha_s - hyperparameter for positive rates theta_l,...

In oslerinhealth/rewind: Reconstructing Etiology with Binary Decomposition

Description

Function to sample the parameters from a grid. NB: can just fix it at a particular value.

Usage

update_pseudo_n_PR(
  theta_vec,
  a_theta_one,
  e = 1,
  f = 0.1,
  show_density = FALSE
)

Arguments

theta_vec	vector of positive rates
a_theta_one	scalar between 0 and 1.
e, f	hyperparameter for Gamma distribution over original alpha_s (mean e/f, variance e/f^2).
show_density	show the full conditional density of alpha_s given other unknown parameters; default is FALSE.

Value

an updated alpha_s value (positive)

Examples

s  <- 10#rgamma(1,0.1,0.01)
a1 <- 0.9
n  <- 20
p  <- rbeta(n, s*a1,s*(1-a1))

par(mfrow=c(1,3))
plot(p,main="raw data: probabilities")

update_pseudo_n_PR(p,a1,e = 0.1,f=0.01,show_density =TRUE)
abline(v=s/(1+s),col="blue")
legend("topleft",c("True s/(1+s)","simulated s/(1+s)"),
       col=c("blue","red"),lty=c(1,2))
mtext(paste0("true s: ",as.character(s)),3)

update_mu_PR(p,s,a=2,b=1,show_density=TRUE)
abline(v=a1,col="blue")
legend("topleft",c("True a_theta1","simulated a_theta1"),
       col=c("blue","red"),lty=c(1,2))
mtext(paste0("true a_theta1: ",as.character(a1)),3)

## Not run: 
par(mfrow=c(1,2))
a10 <- 0.9
s0  <- 5
miter <- 100
res <- matrix(NA,nrow=2,ncol=miter)
res[,1] <- c(a10,s0)
for (iter in 2:miter){
  res[2,iter] <- update_pseudo_n_PR(p,res[1,iter-1],e = 0.1,f=0.01,show_density =!TRUE)
  res[1,iter] <- update_mu_PR(p,res[2,iter],a=2,b=1,show_density=!TRUE)
}
plot(res[1,],type="l",main="mu")
abline(h=c(a1),col=c(2))
plot(res[2,],type="l",main="s")
abline(h=c(s),col=c(2))

## End(Not run)

oslerinhealth/rewind documentation built on May 26, 2021, 6:56 a.m.