Biased coin, unknown number of tosses

Based on an example by Ben Lambert in “A student’s guide to Bayesian statistics”.

1 Description of experiment

A biased coin is being tossed (probability of heads \(\theta > 0.5\)). The coin is tossed \(n\) times, recording the number of heads \(X_i\). The entire experiment is repeated 50 times. In every try the parameters \(\theta\) and \(n\) are the same.

2 Model definition:

\[ X_i \sim binomial(n, \theta) \] \[ n \sim \text{discrete-uniform}(11,20) \] \[ \theta \sim beta(7,2) \] Priors:

• Assuming the number of tosses can be any integer number between 11 and 20.
  
• Assuming the coin is biased for higher number of heads, the probability of heads is ~0.78.

Binomial and Beta distributions:

par(mfrow=c(1,3))
plot(dbinom(c(0:13),13,0.6),type="h",xlab="x",ylab="density",main="binomial(13,0.6)")
plot(dunif(c(1:25),11,020),type="h",xlab="x",ylab="density",main="uniform(11,20)")

xticks<-seq(0,1,length=50)
plot(x=xticks, y=dbeta(xticks,7,2), type="l",xlab="x",ylab="density",main="beta(7,2)")

3 The code

Prepare sample data and load the model. The model will only see the list of numbers \(X\)

library(rstan)
set.seed(2020)
X<-rbinom(50, 13, 0.6)
K<-length(X)
shift<-10
aModel<-stan_model('~\\R\\bayes-talk\\coinFlip.stan')
X

 [1]  7  8  7  8 10 10 10  8 12  7  7  7  6  8  8  8  5  7  8  9  9 10  6  5  6 10  8  7  8  8  5 10  5  5 13  8 12  5  8  9  4  9  6  6
[45]  7  8  8  7  6 10

Run sampler

options(mc.cores=4)
fit<-sampling(aModel,data=list(X=X,K=K),iter=1000,chains=4)

Extract model parameters, estimate \(\theta\)

fit_states<-subset(dimnames(fit)$parameters, grepl("^pState",dimnames(fit)$parameters))
fit_ss<-extract(fit, pars=c('theta',fit_states))
lmean<-lapply(fit_ss, mean)
(mtheta<-lmean$theta)

[1] 0.5082172

quantile(fit_ss$theta,probs=c(0.05,0.95))

       5%       95% 
0.4030376 0.6067676 

Estimate the number of tosses \(n\)

maxpstate<-which.max(as.numeric(lmean[fit_states]))
(maxpstate<-maxpstate+shift)

[1] 15

Plot model parameters space

par(mfrow=c(2,5))
ylim<-c(0,1)
for (i in 1:length(fit_states)) {
  with(fit_ss, {
    plot(x=theta, ylim=ylim, ylab=paste0("p(Tosses=",i+shift,")"),
         main=paste0("p=",round(as.numeric(lmean[fit_states[i]]),digits=2)), 
         pch=".", y=fit_ss[[fit_states[i]]])
    abline(v=quantile(theta,probs=c(0.05,0.5,0.95)), col="gray")
  })
}

4 Analysis

Properties of binomial distribution:

• Probability that a count of heads equals x, given a number of tosses n and a probability of heads p: \[ prob(x) = \binom{n}{x}p^x(1-p)^{n-x} \]

• Expected value: \[ \mathrm{E}(X)=np \]

• Standard deviation: \[ sd(X)=\sqrt{np(1-p)} \]

Given these, let’s evaluate model fit by comparing model-predicted values of sample mean and standard deviation against what the sample has and what the theory provides

maxpstate*mtheta

[1] 7.623259

mean(X)

[1] 7.76

13*0.6

[1] 7.8

sqrt(maxpstate*mtheta*(1-mtheta))

[1] 1.93623

sd(X)

[1] 1.964584

sqrt(13*0.6*0.4)

[1] 1.766352

Differences are explained by small sample size. On larger samples all estimated parameters converge to true model parameters.