R2jags
(meetup pizza example)My Journey From Frequentist to Bayesian Statistics
When there is WinBUGS already, why is JAGS necessary?
install.packages("R2jags")library(R2jags)
Changing our beliefs based on evidences Posterior=Likelihood×PriorEvidence
Nat Napoletano:Bayes' Theorem for Everyone
The focus of the analysis presented here is on accurately estimating the mean of IQ using simulated data.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4158865/
https://towardsdatascience.com/visualizing-beta-distribution-7391c18031f1
We hold meetings every month and post it on meetup.com, which tells us how many people've signed up, n. But not everyone show up (i.e., traffic). So, we'd like to know what percentage of people, θ, will show up to order enough pizzas.
Yi|θ∼Binom(ni,θ)
θ∼Beta(a,b)
Weak prior:
Beta(1,1)↔Unif(0,1)
Y <- c(77, 62, 30)N <- c(39, 29, 16)
But as you know, someone (@lukanegoita) said
Yi|θ∼Binom(ni,θ)
θ∼Beta(a,b)
In R
model{for(i in 1:3){ Y[i] <- rbinom(1, N[i], theta) } theta <- rbeta(1, 1, 1)}
In JAGS
model{for(i in 1:3){ Y[i] ~ dbinom(theta, N[i])} theta ~ dbeta(1, 1)}
dist
instead of rdist
for random samplingsuppressMessages(library(R2jags))model_string <-"model{# Likelihoodfor(i in 1:3){ Y[i] ~ dbinom(theta, N[i])}# Prior theta ~ dbeta(1, 1)}"# DataN <- c(77, 62, 30)Y <- c(39, 29, 16)# specify parameters to be moniteredparams <- c("theta")
R2jags
and runjags
are just wrappers of rjagsrunjags
makes it easy to run chains in parallel on multiple processesR2jags
allows running until converged and paralleling MCMC chainsR2jags
output <- jags(model.file = textConnection(model_string), data = list(Y = Y, N = N), parameters.to.save = params, n.chains = 2, n.iter = 2000)
## Compiling model graph## Resolving undeclared variables## Allocating nodes## Graph information:## Observed stochastic nodes: 3## Unobserved stochastic nodes: 1## Total graph size: 9## ## Initializing model## ## | | | 0% | |** | 4% | |**** | 8% | |****** | 12% | |******** | 16% | |********** | 20% | |************ | 24% | |************** | 28% | |**************** | 32% | |****************** | 36% | |******************** | 40% | |********************** | 44% | |************************ | 48% | |************************** | 52% | |**************************** | 56% | |****************************** | 60% | |******************************** | 64% | |********************************** | 68% | |************************************ | 72% | |************************************** | 76% | |**************************************** | 80% | |****************************************** | 84% | |******************************************** | 88% | |********************************************** | 92% | |************************************************ | 96% | |**************************************************| 100%
## Compiling model graph## Resolving undeclared variables## Allocating nodes## Graph information:## Observed stochastic nodes: 3## Unobserved stochastic nodes: 1## Total graph size: 9## ## Initializing model## ## | | | 0% | |** | 4% | |**** | 8% | |****** | 12% | |******** | 16% | |********** | 20% | |************ | 24% | |************** | 28% | |**************** | 32% | |****************** | 36% | |******************** | 40% | |********************** | 44% | |************************ | 48% | |************************** | 52% | |**************************** | 56% | |****************************** | 60% | |******************************** | 64% | |********************************** | 68% | |************************************ | 72% | |************************************** | 76% | |**************************************** | 80% | |****************************************** | 84% | |******************************************** | 88% | |********************************************** | 92% | |************************************************ | 96% | |**************************************************| 100%
plot(as.mcmc(output))
library(mcmcplots)mcmcplot(as.mcmc(output),c("theta"))denplot(as.mcmc(output),c("theta"))
So it might take a longer time for run
output$BUGSoutput$summary
## mean sd 2.5% 25% 50% 75%## deviance 14.6685495 1.42564249 13.6497484 13.7569065 14.104791 15.0374241## theta 0.4962242 0.03867905 0.4226714 0.4693877 0.496706 0.5209497## 97.5% Rhat n.eff## deviance 18.7817874 1.002871 630## theta 0.5723382 1.000970 2000
If it fails to converge, we can do the following,
update(output, n.iter = 1000)autojags(output, n.iter = 1000, Rhat = 1.05, n.update = 2)
output$BUGSoutput$summary
## mean sd 2.5% 25% 50% 75%## deviance 14.6685495 1.42564249 13.6497484 13.7569065 14.104791 15.0374241## theta 0.4962242 0.03867905 0.4226714 0.4693877 0.496706 0.5209497## 97.5% Rhat n.eff## deviance 18.7817874 1.002871 630## theta 0.5723382 1.000970 2000
output2$BUGSoutput$summary
## mean sd 2.5% 25% 50% 75%## deviance 15.2204981 1.7165810 13.6500147 13.9259913 14.6814049 15.9249686## theta 0.5346611 0.0298314 0.4752045 0.5146234 0.5355181 0.5547805## 97.5% Rhat n.eff## deviance 19.846200 1.000531 2000## theta 0.591993 1.000510 2000
denplot(as.mcmc(output), parms = "theta", xlim = c(0.3, 0.7))
denplot(as.mcmc(output2), parms = "theta", xlim = c(0.3, 0.7))
In JAGS
model<-"model{# Likelihoodfor(i in 1:N){ y[i] ~ dnorm(mu[i], tau) mu[i] <- beta0 + beta1*x[i]}# Prior # intercept beta0 ~ dnorm(0, 0.0001) # slope beta1 ~ dnorm(0, 0.0001) # error term tau ~ dgamma(0.01, 0.01) sigma <- 1/sqrt(tau)}"
In statistics 101
Yi=b0+b1Xi+σ Y∼N(b0+b1X,σ2)
In R
lm(y~x)
## Compiling model graph## Resolving undeclared variables## Allocating nodes## Graph information:## Observed stochastic nodes: 15## Unobserved stochastic nodes: 3## Total graph size: 74## ## Initializing model
In JAGS
lm_output$BUGSoutput$summary[1:2, c(3,7)]
## 2.5% 97.5%## beta0 -8.853992 52.131691## beta1 -16.496277 3.192873
95% credible interval
In R
confint(lm(y~x))
## 2.5 % 97.5 %## (Intercept) -9.32782 51.322335## x -16.19329 3.486528
95% confidence interval
RJAGS: how to get started by Rens van de Schoot, retweeted by Martyn Plummer
"perky ears" + "floppy ears" ?= "half-up ears"
For those still can't decide which any external tool?
Keyboard shortcuts
↑, ←, Pg Up, k | Go to previous slide |
↓, →, Pg Dn, Space, j | Go to next slide |
Home | Go to first slide |
End | Go to last slide |
Number + Return | Go to specific slide |
b / m / f | Toggle blackout / mirrored / fullscreen mode |
c | Clone slideshow |
p | Toggle presenter mode |
t | Restart the presentation timer |
?, h | Toggle this help |
Esc | Back to slideshow |