Question One - Sheep Tick Data

i. Maximum Likelihood Estimation for Lambda

We have:
\[f(y\mid\lambda)=\dfrac{e^{-\lambda}\lambda^y}{y!}\]
Writing out the likelihood for this:
\[L(y\mid\lambda)=\dfrac{e^{-\lambda}\lambda^{y_1}}{y_1!}\cdot\dfrac{e^{-\lambda}\lambda^{y_2}}{y_2!}\cdot...\cdot\dfrac{e^{-\lambda}\lambda^{y_n}}{y_n!}\]
Simplifying:
\[L(y\mid\lambda)=\dfrac{e^{-n\lambda}\lambda^{\sum_{i=1}^ny_i}}{\sum_{i=1}^ny_1!}\]
\[L(y\mid\lambda)={e^{-n\lambda}\cdot\lambda^{\sum_{i=1}^ny_i}} \cdot\frac{1}{\sum_{i=1}^ny_i!}\]
Taking the log likelihood:
\[l(y\mid\lambda)=ln({e^{--n\lambda}\cdot\lambda^{\sum_{i=1}^ny_i}} \cdot\frac{1}{\sum_{i=1}^ny_i!})\]
\[l(y\mid\lambda)=(-n\lambda)+ln(\lambda^{\sum_{i=1}^ny_i}) + ln(\frac{1}{\sum_{i=1}^ny_i!})\] \[l(y\mid\lambda)=(-n\lambda)+({\sum_{i=1}^ny_i})ln(\lambda) -ln({\sum_{i=1}^ny_i!})\]
Differentiating with respect to \(\lambda\) \[\dfrac{\partial l}{\partial \lambda}=-n+\dfrac{\sum_{i=1}^ny_i}{\lambda}\]
Setting the derivative to 0 and solving for \(\lambda\)
\[-n+\dfrac{\sum_{i=1}^ny_i}{\lambda}=0\]
\[\dfrac{\sum_{i=1}^ny_i}{\lambda}=n\]
\[\sum_{i=1}^ny_i=n\lambda\]
Since \[\sum_{i=1}^ny_i=n\bar{y}\]
\[n\bar{y}=n\lambda\]
\[\frac{n\bar{y}}{n}=\lambda\]
\[\hat{\lambda}= \bar{y}\]
Now verifying the stationary point was indeed a maximum by taking second derivative:
\[\dfrac{\partial l}{\partial \lambda}=-n+\dfrac{\sum_{i=1}^ny_i}{\lambda}\]
\[\dfrac{\partial^2 l}{\partial \lambda^2}=-\dfrac{\sum_{i=1}^ny_i}{\lambda^2}\leq0\]
As the poisson random variables \(y_i\) and \(\lambda^2\) will always be positive, the negative one scalar guarantees that the stationary point is a minimum. Thus the desired result have been achieved, verifying that:
MLE \(\hat{\lambda}= \bar{y}\)

ii. Finding the Variance of Lambda

\[var({\hat{\lambda}})=var({\bar{y}})\]
\[var({\hat{\lambda}})=var(\frac{\sum^n_{i=1}y_i}{n})\]
\[var({\hat{\lambda}})=\frac{1}{n^2}\cdot var({\sum^n_{i=1}y_i})\]
\[var(Y)=\lambda\]
\[var({\hat{\lambda}})=\frac{1}{n^2}\cdot n\lambda)\]
\[var({\hat{\lambda}})=\frac{\hat{\lambda}}{n}\]

iii. Frequentist Estimation of Lambda From Fisher’s Tick Data

library(kableExtra)
# counts of ticks (col 1 is the number of ticks, col 2 is the number of sheep with that many ticks)
xcount = matrix(c(0, 4,
1, 5,
2, 11,
3, 10,
4, 9,
5, 11,
6, 3,
7, 5,
8, 3,
9, 2,
10, 2,
11, 5,
12, 0,
13, 2,
14, 2,
15, 1,
16, 1,
17, 0,
18, 0,
19, 1,
20, 0,
21, 1,
22, 1,
23, 1,
24, 0,
25, 2),26,2,byrow=T) # 
y = rep(xcount[,1],times=xcount[,2]) # gives a vector x_1,...,x_n of number of ticks per sheep
lambda_hat<-mean(y) #calculating lambda hat using previous MLE
se<-lambda_hat/length(y) #standard error using previous result
z<-qnorm(0.975) #our z statistic
ci_l<-lambda_hat-z*se #lower bound
ci_u<-lambda_hat+z*se #upper bound
ci<-c(ci_l,ci_u)
freq1<-(rbind(c(lambda_hat,""), ci))
rownames(freq1) =c("Lambda_Hat/Mean", "95% Confidence Interval")
kable(freq1, align=c("l" , "c","l"), caption="Frequentist Inference")  |>kable_classic_2(full_width = F, html_font = "Cambria")
Frequentist Inference
Lambda_Hat/Mean 6.5609756097561
95% Confidence Interval 6.40415517196869 6.71779604754351

Our maximum likelihood estimate for lambda is 6.560976. Our 95% confidence interval has bounds of 6.404155 and 6.717796, meaning 95% of the time we expect the true value of lambda to lie between this range.

Using the previous frequentist results with the sheep tick data, \(\hat{\lambda}\) is found to be 6.560976. The 95% confidence interval for \(\lambda\) is (6.404155, 6.717796), meaning that 95% of the time we would expect the true value of \(\lambda\) to lie within this range.

iv. Expressing the Posterior Distribution of Lambda

First we will define our likelihood/data (using what was found earlier for the MLE calculation) and prior distribution:
\[f(y_1,y_2,...,y_n\mid \lambda)=e^{-n\lambda}\cdot\lambda^{\sum^n_{i=1}y_i}\cdot\frac{1}{\sum^n_{i=1}y_i !}\]
\[f(\lambda)=\frac{\beta^{\alpha} }{\Gamma(\alpha)}\cdot \lambda^{\alpha - 1} e^{-\beta \lambda}\]
where:
\(\alpha > 0\) is the shape parameter,
\(\beta > 0\) is the rate parameter,
\(\Gamma(\alpha)\) is the Gamma function.

By Bayes’ theorem we can write (then substitute and simplify):
\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto f(y_1,y_2,...,y_n\mid \lambda) \, f(\lambda)\]

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha)}\cdot \lambda^{\alpha - 1} e^{-\beta \lambda} \cdot e^{-n\lambda}\cdot\lambda^{\sum^n_{i=1}y_i}\cdot\frac{1}{\sum^n_{i=1}y_i !}\]

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha)}\cdot e^{-\beta \lambda-n\lambda} \cdot\frac{\lambda^{\alpha - 1}\cdot \lambda^{\sum^n_{i=1}y_i}}{\sum^n_{i=1}y_i !}\]

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha)}\cdot e^{-(\beta +n)\lambda} \cdot\frac{\lambda^{((\alpha+\sum^n_{i=1}y_i) - 1)}}{\sum^n_{i=1}y_i !}\]

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha) \sum^n_{i=1}y_i !}\cdot e^{-(\beta +n)\lambda} \cdot{\lambda^{((\alpha+\sum^n_{i=1}y_i) - 1)}}\]

Now I will define two new parameters:
\[\boldsymbol{a} = \alpha+\sum^n_{i=1}y_i\]
\[\boldsymbol{b} = \beta +n\]

Further simplifying my posterior:
\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha) \sum^n_{i=1}y_i !}\cdot e^{-\boldsymbol{b}\lambda} \cdot{\lambda^{(\boldsymbol{a} - 1)}}\]

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha) \sum^n_{i=1}y_i !}\cdot {\lambda^{(\boldsymbol{a} - 1)}} \cdot e^{-\boldsymbol{b}\lambda} \]

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \overbrace{\frac{\beta^{\alpha} }{\Gamma(\alpha)} \underbrace{\frac{1}{\sum^n_{i=1}y_i !}}_\text{Constant} \cdot \underbrace{{\lambda^{(\boldsymbol{a} - 1)}} \cdot e^{-\boldsymbol{b}\lambda}}_\text{Kernel of Gamma Dist.} }^\text{Gamma Posterior Distribution} \newline\]
\[\boldsymbol{a}=\alpha+\sum^n_{i=1}y_i\] \[\newline \boldsymbol{b}=\beta +n\]

Full Simplified Posterior:

\[f(\lambda \mid {y_1,y_2,...,y_n}) \propto \frac{\beta^{\alpha} }{\Gamma(\alpha)}\cdot {\lambda^{(\boldsymbol{a} - 1)}} \cdot e^{-\boldsymbol{b}\lambda} \]

v. Bayesian Estimation for the value Lambda

alpha<-0.001
beta<-0.001
a_1<-alpha+sum(y)
b_1<-beta+length(y)
mean_post<-a_1/b_1
median_post<-qgamma(0.5,a_1,b_1)
credible_interval<-qgamma(c(0.025,0.975),a_1,b_1)
library(knitr)
library(kableExtra)
cool_stuff<-rbind("Posterior Mean"=c(mean_post, ""),"Posterior Median"= c(median_post, ""),"MLE Lambda"=c(lambda_hat, ""), "95% Credible Interval" =credible_interval)
kable(cool_stuff, align=c("l" , "c"), caption="Sanity Check Statistics")|>kable_classic_2(full_width = F, html_font = "Cambria")
Sanity Check Statistics
Posterior Mean 6.56090779380739
Posterior Median 6.55684325078098
MLE Lambda 6.5609756097561
95% Credible Interval 6.01815062613793 7.12676288167517

From our posterior distribution, we find the mean value of Lambda to be 6.560908, while the median estimate is 6.556843. This suggest a centred distribution with a very slight right skew, characteristic of our gamma posterior given the shape and rate parameters. It is also reassuring to find these estimate close to our maximum likelihood result of \(\hat{\lambda}\) = 6.560976. Our 95% credible interval ranges between 6.018151 and 6.126763, meaning there is a 95% probability that \(\hat{\lambda}\) falls within this range based on our posterior distribution.

vi. JAGS Model

Markov Chains

library(rjags)
## Loading required package: coda
## Linked to JAGS 4.3.2
## Loaded modules: basemod,bugs
library(bayesplot)
## This is bayesplot version 1.13.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(posterior)
## This is posterior version 1.6.1
## 
## Attaching package: 'posterior'
## The following object is masked from 'package:bayesplot':
## 
##     rhat
## The following objects are masked from 'package:stats':
## 
##     mad, sd, var
## The following objects are masked from 'package:base':
## 
##     %in%, match
library(ggplot2)

tmp="
model{
  for(i in 1:n){
    y[i]~dpois(lambda)}
  lambda~dgamma(0.001,0.001)
}
"
cat(tmp, file = "model_pois.txt")

n<-as.data.frame(length(y))
init_fun_norm1 = function(){
init_lambda = rnorm(1, mean = 0.5, sd = 0.1)
return(list(lambda = init_lambda))}
m1<-jags.model(file="model_pois.txt", data=list(y=y, n=length(y)), inits = init_fun_norm1, n.chains = 3, n.adapt = 1000)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 82
##    Unobserved stochastic nodes: 1
##    Total graph size: 85
## 
## Initializing model
out1 = coda.samples(m1,c('lambda'),10000)
str(out1)
## List of 3
##  $ : 'mcmc' num [1:10000, 1] 6.62 6.29 6.72 6.92 6.56 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "lambda"
##   ..- attr(*, "mcpar")= num [1:3] 1 10000 1
##  $ : 'mcmc' num [1:10000, 1] 6.1 6.71 7.08 6.45 6.17 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "lambda"
##   ..- attr(*, "mcpar")= num [1:3] 1 10000 1
##  $ : 'mcmc' num [1:10000, 1] 6.64 6.89 6.6 6.58 6.43 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "lambda"
##   ..- attr(*, "mcpar")= num [1:3] 1 10000 1
##  - attr(*, "class")= chr "mcmc.list"
out1_df = as_draws_df(out1)
out1_df

Posterior Distribution Graph

ggplot(out1_df, aes(x=lambda))+geom_histogram(aes(y=after_stat(density),fill="Histogram"),color="pink", bins=50)+labs(title="Histogram of Posterior Density",x="Lambda Value", y = "Density")+ geom_vline(aes(xintercept=mean(lambda), colour="Posterior Mean"),linetype= "dashed")+ stat_function(fun = dgamma, args = list(shape = a_1, rate = b_1),aes( color = "Gamma Density"), linetype = "dashed") +  scale_fill_manual(name = "", values = c("Histogram" = "powderblue")) +scale_color_manual(name = "", values = c("Posterior Mean" = "black","Gamma Density" = "purple")) +
theme_minimal() +guides(fill = guide_legend(order = 1), 
         color = guide_legend(order = 2))

Corroborating with the results from our manual calculation of the posterior distribution, we have a nearly symmetrical distribution centered at a mean of around 6.55. There is a slight right skew. The blue represents the density of 10,000 samples of Lambda from our posterior distribution and given parameters, while the purple represents gamma distribution with the same posterior parameters.

Markov Chains Overlay Graph

mcmc_dens_overlay(out1_df, pars = "lambda")

A distribution showing the density of each of the 3 Markov chains used in forming our posterior.

Table of Results

credible_interval<-quantile(out1_df$lambda,c(0.025,0.975))
c<-kable(credible_interval, col.names=c("Quantile", "Lambda Value"), align=c("l" , "c","l"), caption="95% Credible Interval for Lambda")  |>kable_classic_2(full_width = F, html_font = "Cambria")
c
95% Credible Interval for Lambda
Quantile Lambda Value
2.5% 6.021557
97.5% 7.133773 
sout1_r<-as.data.frame(t(summary(out1_df))) 
rownames(sout1_r)=c("Variable", "Mean", "Median","Standard Deviation", "Mean Absolute Deviation", "Lower 5% Quantile", "Upper 5% Quantile", "R-Hat", "Ess_Bulk", "Ess_Tail")
kable(sout1_r, col.names = c("Statistic", "Value"), align=c("l" , "c","l"), caption="Summary Statistics from JAGS")  |>kable_classic_2(full_width = F, html_font = "Cambria")
Summary Statistics from JAGS
Statistic Value
Variable lambda
Mean 6.562894
Median 6.559693
Standard Deviation 0.2830349
Mean Absolute Deviation 0.282614
Lower 5% Quantile 6.106625
Upper 5% Quantile 7.037798
R-Hat 1.000004
Ess_Bulk 30207.19
Ess_Tail 29228.74

From the JAGS output we have a 95% credible interval for the value of lambda with a lower bound of 6.0215568 and upper bound of 7.1337731 . The probability that the true value of Lambda falls within this range in 95%.

The mean is 6.562894 and the median 6.559693, both very similar to what was found with the manual computation of the posterior distribution. The standard deviation of lambda is 0.2830349. The R-Hat value is very close to 1 with 1.000004, indicating convergence of the chains.

Trace Plot

color_scheme_set("viridis")
mcmc_trace(out1_df, pars = "lambda")

Looking at the trace plot we see significant overlap between the chains with close to no drifting. This is ‘good behaviour’, giving confidence the chains have reached their limiting distribution. The convergence is corroborated by our R-Hat value.

Question 2 - Gull Data

i. Finding the Posterior Distrubtion for Pi

First we will find and define the likelihood/data, followed by the prior distribution: \[f(y_1,y_2,...,y_n \mid \pi)=\frac{n!}{y_1!(n-y_1)!}\pi^{y_1}(1-\pi)^{n-y_1}\cdot\frac{n!}{y_2!(n-y_2)!}\pi^{y_2}(1-\pi)^{n-y_2}\cdot...\cdot\frac{n!}{y_n!(n-y_n)!}\pi^{y_n}(1-\pi)^{n-y_n}\]
\[f(y_1,y_2,...,y_n \mid \pi)=\prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\pi^{y_i}(1-\pi)^{n-y_i}\] \[f(y_1,y_2,...,y_n \mid \pi)=\prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\cdot\pi^{\sum_{i=1}^ny_i}\cdot(1-\pi)^{n-\sum_{i=1}^ny_i}\]

\[f(\pi)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \pi^{\alpha - 1} (1-\pi)^{\beta-1}\]

where:
\(\alpha > 0\) is the shape parameter,
\(\beta > 0\) is the rate parameter,
\(\Gamma(\alpha)\) is the Gamma function.

By Bayes’ theorem we can write (then substitute and simplify): \[f(\pi \mid {y_1,y_2,...,y_n}) \propto f(y_1,y_2,...,y_n\mid pi) \, f(\pi)\]
\[f(\pi \mid {y_1,y_2,...,y_n}) \propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \pi^{\alpha - 1} (1-\pi)^{\beta-1}\cdot \prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\cdot\pi^{\sum_{i=1}^ny_i}\cdot(1-\pi)^{n-\sum_{i=1}^ny_i}\]

\[f(\pi \mid {y_1,y_2,...,y_n}) \propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\cdot \pi^{\alpha - 1} (1-\pi)^{\beta-1}\cdot \pi^{\sum_{i=1}^ny_i}\cdot(1-\pi)^{n-\sum_{i=1}^ny_i}\]
\[f(\pi \mid {y_1,y_2,...,y_n}) \propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\cdot \pi^{\alpha - 1}\cdot \pi^{\sum_{i=1}^ny_i} (1-\pi)^{\beta-1}\cdot(1-\pi)^{n-\sum_{i=1}^ny_i}\] \[f(\pi \mid {y_1,y_2,...,y_n}) \propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\cdot \pi^{\alpha - 1+\sum_{i=1}^ny_i}\cdot(1-\pi)^{\beta-1+n-\sum_{i=1}^ny_i}\] Now I will define two new parameters:
\[\boldsymbol{a} = \alpha+\sum^n_{i=1}y_i\]
\[\boldsymbol{b} = \beta +n - \sum^n_{i=1}y_i\]
\[f(\pi \mid {y_1,y_2,...,y_n}) \propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}\cdot \pi^{\boldsymbol{a} - 1}\cdot(1-\pi)^{\boldsymbol{b}-1}\]
\[f(\pi \mid {y_1,y_2,...,y_n}) \propto \overbrace{\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \underbrace{\cdot \prod_{i=1}^n\frac{n!}{y_i!(n-y_i)!}}_\text{Constant} \cdot \underbrace{\cdot \pi^{\boldsymbol{a} - 1}\cdot(1-\pi)^{\boldsymbol{b}-1} \cdot e^{-\boldsymbol{b}\lambda}}_\text{Kernel of Beta Dist.}}^\text{Beta Posterior Distribution} \newline\] \[\boldsymbol{a} = \alpha+\sum^n_{i=1}y_i\]
\[\boldsymbol{b} = \beta +n - \sum^n_{i=1}y_i\] Fully simplified posterior:
\[f(\pi \mid {y_1,y_2,...,y_n}) \propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot \pi^{\boldsymbol{a} - 1}\cdot(1-\pi)^{\boldsymbol{b}-1}\]

ii. Bayesian Estimation of the Value of Pi

n_2<-609
d<-175
alpha<-0.5
beta<-0.5
a_2<-alpha+d
b_2<-beta+n_2-d
mean_post2<-a_2/(a_2+b_2)
median_post2<-qbeta(0.5,a_2,b_2)
credible_intervalB<-qbeta(c(0.025,0.975),a_2,b_2)
library(knitr)
library(kableExtra)
cool_stuff2<-rbind("Posterior Mean"=c(mean_post2, ""),"Posterior Median"= c(median_post2, ""), "95% Credible Interval" =credible_intervalB)
kable(cool_stuff2, align=c("l" , "c"), caption="Sanity Check Statistics")|>kable_classic_2(full_width = F, html_font = "Cambria")
Sanity Check Statistics
Posterior Mean 0.287704918032787
Posterior Median 0.287472769434021
95% Credible Interval 0.252486363488668 0.324241989509822

iii. Probability of Pi Lower than 0.3

pi_03<-pbeta(0.3,a_2, b_2, lower.tail=TRUE )
kable(pi_03, align=c("l" , "c","l"), col.names = ("Probability"), caption="Probability of Pi Less Than 0.3")  |>kable_classic_2(full_width = F, html_font = "Cambria")
Probability of Pi Less Than 0.3
Probability
0.7510092

The probability of observing a value of 0.3 or lower for pi is 0.7510092.

iv. JAGS Model

Markov Chains

tmp3="
model{
    d~dbin(pi,n)
  pi~dbeta(0.5,0.5)
}
"
cat(tmp3, file = "model_binom.txt")

n<-n_2
init_fun_norm2 = function(){
init_pi = rnorm(1, mean = 0.5, sd = 0.1)
return(list(pi = init_pi))}
m_b<-jags.model(file="model_binom.txt", data=list(d=d, n=n_2), inits = init_fun_norm2, n.chains = 3, n.adapt = 1000)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 4
## 
## Initializing model
outb = coda.samples(m_b,c('pi'),10000)
str(outb)
## List of 3
##  $ : 'mcmc' num [1:10000, 1] 0.33 0.299 0.287 0.273 0.284 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pi"
##   ..- attr(*, "mcpar")= num [1:3] 1001 11000 1
##  $ : 'mcmc' num [1:10000, 1] 0.298 0.274 0.296 0.28 0.267 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pi"
##   ..- attr(*, "mcpar")= num [1:3] 1001 11000 1
##  $ : 'mcmc' num [1:10000, 1] 0.294 0.263 0.299 0.286 0.316 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pi"
##   ..- attr(*, "mcpar")= num [1:3] 1001 11000 1
##  - attr(*, "class")= chr "mcmc.list"
outb_df = as_draws_df(outb)
outb_df

Histogram of the Posterior Distribution for Pi

ggplot(outb_df, aes(x=pi))+geom_histogram(aes(y=after_stat(density),fill="Histogram"),color="brown", alpha=0.7, bins= 50)+labs(title="Histogram of Posterior Density",x="Pi Value", y = "Density")+ geom_vline(aes(xintercept=mean(pi), colour="Posterior Mean"),linetype= "dashed")+ stat_function(fun = dbeta, args = list(shape = a_2, shape2 = b_2),aes( color = "Beta Density"), linetype = "dashed") +  scale_fill_manual(name = "", values = c("Histogram" = "cornflowerblue")) +scale_color_manual(name = "", values = c("Posterior Mean" = "red","Beta Density" = "black")) +
theme_minimal() +guides(fill = guide_legend(order = 1), 
         color = guide_legend(order = 2))

The histogram formed from 10,000 posterior distribution samples of pi is roughly symmetrical and bell shaped, corroborating with a beta posterior.. The mean is around 0.285 while the median appears closer to 0.28. Superimposed is the true beta density with the posterior parameters calculated earlier, agreeing with our density of pi samples.

Markov Chains Overlay

mcmc_dens_overlay(outb_df, pars = "pi")

An overlay of each Markov chain.

Table of Results

credible_interval_mb<-quantile(outb_df$pi,c(0.025,0.975))
find_03<-ecdf(outb_df$pi)(0.3)
kable(credible_interval_mb, col.names=c("Quantile", "Pi Value"), align=c("l" , "c","l"), caption="95% Credible Interval for Pi")  |>kable_classic_2(full_width = F, html_font = "Cambria")
95% Credible Interval for Pi
Quantile Pi Value
2.5% 0.2526947
97.5% 0.3236652 
sout2_r<-rbind(as.data.frame(t(summary(outb_df))), find_03) 
rownames(sout2_r)=c("Variable", "Mean", "Median","Standard Deviation", "Mean Absolute Deviation", "Lower 5% Quantile", "Upper 5% Quantile", "R-Hat", "Ess_Bulk", "Ess_Tail", "P(pi<0.3)")
kable(sout2_r, col.names = c("Statistic", "Value"), align=c("l" , "c","l"), caption="Summary Statistics from JAGS")  |>kable_classic_2(full_width = F, html_font = "Cambria")
Summary Statistics from JAGS
Statistic Value
Variable pi
Mean 0.2875359
Median 0.2872026
Standard Deviation 0.01818988
Mean Absolute Deviation 0.01817755
Lower 5% Quantile 0.2578973
Upper 5% Quantile 0.3178528
R-Hat 1.000022
Ess_Bulk 18841.19
Ess_Tail 17557.22
P(pi<0.3) 0.754666666666667

The 95% credible interval for pi has a lower bound of 0.2526947 and upper bound of 0.3236652 . The probability of the true value of pi occurring between these two bounds is 95%. Looking at the central tendency of our posterior distribution, the mean value of pi is 0.2875359 while the median is marginally lower at 0.2872026. The standard deviation is 0.01818988. The R-hat value is very close to 1 at 1.000022, giving confidence that the Markov Chains converge. Now using the JAGS output, the probability of observing a pi value less than 0.3, meaning 30% or less of the gulls divorce, is `0.75190.7546667 which is marginally higher than the manual calculation.

Overall this density found with JAGS is very similar to the true density found earlier. The 95% credible interval, mean and median are nearly identical.

Trace Plot

color_scheme_set("viridis")
mcmc_trace(outb_df, pars = "pi")

The trace plot indicates ‘good behaviour’ of the Markov Chains. There is no drifting, while the chains have significant overlap. Thus the ‘fuzzy catepillar’ appearance indicates the chains have reached their limiting distribution. This convergence is indicated by the R-Hat value close to 1.

Question 3 - Dunedin Weather Data

i. JAGS Model

Markov Chains

y2 <-
c(0, 0, 0, 1, 0, 0, 0, 2, 5, 0, 0, 0, 0, 3, 1, 0, 0, 10, 0, 0,
3, 14, 4, 2, 4, 8, 0, 1, 1, 1, 1, 7, 2, 0, 0, 0, 5, 1, 0, 0,
0, 2, 3, 1, 0, 0, 0, 1, 0, 3, 1, 0, 1, 1, 0, 3, 0, 0, 0, 0, 3,
4, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0)
tmp2="
model{
  for(i in 1:n){
    y[i]~dnegbin(pi,1)}
  pi~dbeta(1,1)
  mean<-1/pi
}
"
cat(tmp2, file = "model_geom.txt")

n2<-as.data.frame(length(y2))
init_fun_norm3 = function(){
init_pi = rnorm(1, mean = 0.5, sd = 0.1)
return(list(pi = init_pi))
} 
m2<-jags.model(file="model_geom.txt", data=list(y=y2, n=length(y2)), inits = init_fun_norm3, n.chains = 3, n.adapt = 1000)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 78
##    Unobserved stochastic nodes: 1
##    Total graph size: 82
## 
## Initializing model
out2 = coda.samples(m2,c('pi','mean'),10000)
str(out2)
## List of 3
##  $ : 'mcmc' num [1:10000, 1:2] 2.83 2.41 2.59 2.46 2.31 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "mean" "pi"
##   ..- attr(*, "mcpar")= num [1:3] 1 10000 1
##  $ : 'mcmc' num [1:10000, 1:2] 2.31 2.77 2.46 2.7 2.02 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "mean" "pi"
##   ..- attr(*, "mcpar")= num [1:3] 1 10000 1
##  $ : 'mcmc' num [1:10000, 1:2] 2.03 2.55 2.27 2.47 2.52 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "mean" "pi"
##   ..- attr(*, "mcpar")= num [1:3] 1 10000 1
##  - attr(*, "class")= chr "mcmc.list"
out2_df = as_draws_df(out2)
out2_df

Histogram of the Posterior Density for Pi

pi_mean <- mean(out2_df$pi)
pi_var <- var(out2_df$pi)
shape1 <- pi_mean * ( (pi_mean * (1 - pi_mean) / pi_var) - 1 )
shape2 <- (1 - pi_mean) * ( (pi_mean * (1 - pi_mean) / pi_var) - 1 )
ggplot(out2_df, aes(x=pi))+geom_histogram(aes(y=after_stat(density),fill="Histogram"),color="darkgreen", bins=50)+labs(title="Histogram of Posterior Density",x="Pi Value", y = "Density")+ geom_vline(aes(xintercept=mean(pi), colour="Posterior Mean"),linetype= "dashed")+ stat_function(fun = dbeta, args = list(shape = shape1, shape2 = shape2),aes( color = "Beta Density"), linetype = "dashed") +  scale_fill_manual(name = "", values = c("Histogram" = "green")) +scale_color_manual(name = "", values = c("Posterior Mean" = "black","Beta Density" = "red")) +
theme_minimal() +guides(fill = guide_legend(order = 1), 
         color = guide_legend(order = 2))

The 10,000 samples of pi from our JAGS model creates this histogram, with densities for different pi values. The mean appears to be approximately 4.4, while the median seems slightly lower at about 4.3. The fitted beta probability density, the red dotted line, indicates how the sample densities follow the posterior distribution.

Markov Chains Overlay

mcmc_dens_overlay(out2_df, pars = "pi")

A density overlay for the value of pi in each chain.

Table of Results

credible_interval2<-quantile(out2_df$pi,c(0.025,0.975))
kable(credible_interval2, col.names=c("Quantile", "Pi Value"), align=c("l" , "c","l"), caption="95% Credible Interval for Pi")  |>kable_classic_2(full_width = F, html_font = "Cambria")
95% Credible Interval for Pi
Quantile Pi Value
2.5% 0.3613722
97.5% 0.5035133 
sout3_r<-as.data.frame(t(summary(out2_df))) 
rownames(sout3_r)=c("Variable", "Mean", "Median","Standard Deviation", "Mean Absolute Deviation", "Lower 5% Quantile", "Upper 5% Quantile", "R-Hat", "Ess_Bulk", "Ess_Tail")
kable(sout3_r, col.names = c("Statistic", "Value"), align=c("l" , "c","l"), caption="Summary Statistics from JAGS")  |>kable_classic_2(full_width = F, html_font = "Cambria")
Summary Statistics from JAGS
Statistic Value
Variable mean pi
Mean 2.3336293 0.4316044
Median 2.3173317 0.4315308
Standard Deviation 0.19988992 0.03630766
Mean Absolute Deviation 0.19513687 0.03643342
Lower 5% Quantile 2.0331055 0.3725467
Upper 5% Quantile 2.6842273 0.4918584
R-Hat 1.000106 1.000120
Ess_Bulk 29614.23 29614.23
Ess_Tail 29243.42 29243.42
The true value of pi has a 95% probability of occurring between a value 0 0.3613722 and 0.5035133, based on the credible interval found. The mean value of pi from the posterior distribution is
, while the median is slightly lower at
. The standard deviation is
. The R-Hat value is very close to 1, at

, which indicates the chains do converge and are reliable.

Trace plot

color_scheme_set("viridis")
mcmc_trace(out2_df, pars = "pi")

The trace plot is a classic ‘fuzzy catepillar’, indicating ‘good behaviour of the chains’. There is no drifting and significant overlap of the chains, indicating the chains have reached their limiting distribution. This is to be expected based on the R-Hat value being so close to 1.

ii. Mean Number of Consecutive Dry Days

Histogram

ggplot(out2_df, aes(x=mean))+geom_histogram(aes(y=after_stat(density)),fill="chocolate",color="darkgreen", bins=50)+labs(title="Histogram of Posterior Density",x="Pi Value", y = "Density")+ geom_vline(aes(xintercept=mean(mean)), colour="black",linetype= "dashed")

The JAGS samples of mean dry days before a rainy day in Dunedin form a somewhat symmetrical distribution with a slight right skew. This is because the mean is equal to \(\frac{1}{\pi}\), thus we have a beta distribution with differing parameters (reminiscent of our previous pi density). The mean appears to occur around 2.3 days while the median is slightly lower.

90% Credible Interval

credible_interval4<-quantile(out2_df$mean,c(0.05,0.95))
kable(credible_interval4, col.names=c("Quantile", "Mean Value"), align=c("l" , "c","l"), caption="90% Credible Interval for the Mean number of Consecutive Dry Days")  |>kable_classic_2(full_width = F, html_font = "Cambria")
90% Credible Interval for the Mean number of Consecutive Dry Days
Quantile Mean Value
5% 2.033106
95% 2.684227

The 90% credible interval for the mean number of consecutive dry days until a rainy day in Dunedin has a lower bound of 2.0331055 and an upper bound of 2.6842273. Thus the true value of the mean has a 90% probability of existing between these bounds. The mean number of these means from our JAGS samples is 2.3336293 while the median is 2.3173317 dry days before a rainy day.

```