Contents

Selected Model

  • Uniform distribution: Discrete and Continuous Uniform distribution
  • Bernoulli process related distribution: Bernoulli, Binomial, Hypergeometric, Negative binomial, Geometric distribution
  • Poisson process related distribution: Poisson (discrete), Exponential, Gamma, Beta
  • Normal distribution related: most of extra models below (Chapter 5)
    • Extra Models (for advanced students’ self study - Those will be included in quiz 3)
      • Chi-squared, T, F, Extreme value, Log-nomal, Cauchy, Weibull, Multinomial, Multivariate normal, Dirichlet distribution

Special references for this chapter

  • http://www.distributome.org/
  • Hand-book on STATISTICAL DISTRIBUTIONS for experimentalists, Walck (4PDFtheoryFull.pdf)
  • Handbook of statistical distribution with applications, Krishnamoorthy, 2006, Chapman & Hall/CRC
  • Mathematical Statistics by Wackerly et al or any others.

Tutorial

Discrete Uniform Distribution (s5)

  • tossing a fair die
# install.packages("prob")
suppressMessages(library(prob))
fairdie <- rolldie(1, makespace = TRUE)
Ex <- sum(fairdie$X1*fairdie$probs)
Ex2 <- sum(fairdie$X1^2*fairdie$probs)
Var_x <- Ex2 - Ex^2
Var_x
## [1] 2.916667
pmf <- fairdie$probs
cdf <- cumsum(pmf)
cdf <- c(0,cdf)

# Plot the probability mass function (pmf)
plot(1:6, fairdie$probs, type = 'h', col = 3, ylim = c(0,0.2), xlab = 'x', ylab = 'f(x)')
points(1:6, fairdie$probs, col = 2)
abline(h = 0, col = 3)

# Plot the cumulative distribution function (cdf)
plot(1:6, cdf[2:length(cdf)], type = 'p', col = 'steelblue', ylim = c(0,1), 
     pch = 19, xlab = 'x', ylab = 'F(x)')
segments(x0 = 0:6, x1 = 0:6 + 1, y0 = cdf)
segments(x0 = 0:6 + 1, y0 = c(cdf[-1], 1), y1 = cdf, lty = 2)
abline(h = 0:1, col = 4)

# alternative approach
# sampling
rdu <- function(n,k) sample(1:k,n,replace = T)

table(rdu(100,6))/100
## 
##    1    2    3    4    5    6 
## 0.25 0.15 0.19 0.10 0.18 0.13
table(rdu(1000,6))/1000
## 
##     1     2     3     4     5     6 
## 0.183 0.140 0.171 0.181 0.173 0.152
table(rdu(10000,6))/10000
## 
##      1      2      3      4      5      6 
## 0.1667 0.1698 0.1669 0.1587 0.1666 0.1713
# pmf
pmf <- function(x,k) ifelse(x >= 1 & x <= k & round(x) == x,1/k,0) 
pmf(1:6,6)
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
# cdf
cdf <- function(x,k) ifelse(x < 1,0,ifelse(x <= k,floor(x)/k,1))
cdf(1:6, 6)
## [1] 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 1.0000000

Continuous uniform distribution (s7)

Suppose that each oyster contains a pearl with a diameter in mm that has a \(U(0,10)\) distribution. Get the expected pearl diamter, variance standard deviation, the probability that the diameter is larger than 3 mm.

a <- 0
b <- 10
x <- 3

# the expected pearl diameter and variance, standard deviation
Ex <- (a + b) / 2
Var <- (b - a)^2 / 12
SD <- sqrt(Var)
Ex; Var; SD
## [1] 5
## [1] 8.333333
## [1] 2.886751
# the probability that the diameter is larger than 3mm
1 - punif(x, min = a, max = b)
## [1] 0.7
x <- seq(-5, 15, length = 1000)
hx <- dunif(x, min = a, max = b)
plot(x, hx, type = "l", lty = 2, xlab = "x value", ylim = c(0,0.2),
     ylab = "Density", main = "uniform distribution")
abline(v = 3, lty = 3, col = "red")

Continuous uniform distribution - Semiconservative replication (s7)

Suppose semiconservative replication occurs every half hour in an organism. What is the probability that the scientist will have to wait at least 20 minutes to observe the semiconservative replication? assume semiconservative replication follows uniform distribution \(U(0,30)\). In other words, the semiconservative replication could equally occur any time between 0 to 30 min. calculate probability of the event occur after 20 min.

a <- 0
b <- 30

# P(X>=20) = 1-P(X<20)
1 - punif(20, min = a, max = b)
## [1] 0.3333333
x <- seq(-5, 35, length = 1000)
hx <- dunif(x, min = a, max = b)
plot(x, hx, type = "l", lty = 2, xlab = "x value", ylim = c(0,0.1),
     ylab = "Density", main = "uniform distribution")
abline(v = 20, lty = 3, col = "red")

Binomial distribution (s11)

Test of toothpaste preventing tooth cavity Incident of cavity in population is 10%. One out of 20 children is found to have new cavities after one year of using new toothpaste.What is the probability to have this result simply by chance? The probability distribution finding children having cavities in general population:

  • \(B(n=20,p=0.1)\)

  • \(P(X=1)\)

n <- 20
p <- 0.1

# sampling 1 child
dbinom(1,n,p)
## [1] 0.2701703
# The test result can happen with ~ 40% of probability by chance. 

plot(0:20, dbinom(0:20, n, p), type = 'h', col = 3, ylim = c(0,0.4), 
     xlab = 'x', ylab = 'f(x)', main = "binomial distribution B(20,0.1)")
points(0:20, dbinom(0:20, n, p), col = 2)
abline(h = 0, col = 3)

Hypergeometric distribution (s13)

If there are 6 defective items out of total 16 items in a box, when 5 items are chosen at random for the inspection, what is the expected number of defective items, and what is the probability to choose two defective items?

N <- 16 # n = N-M
M <- 6 # M
n <- 5 # k
x <- 2 # x

# p(X=2)
dhyper(x, M, N - M, n)
## [1] 0.4120879
# Ex
Ex <- n*M/N
Ex
## [1] 1.875

Hypergeometric distribution - fish in a lake (s13)

An investigator tags 10 fish in a lake containing total 50 fish, and later, catches 8 fish and checks how many fish are tagged, without releasing them until 8 fish are all caught and checked. What is the expected number of tagged fish in this test?

N <- 50 # n = N-M
M <- 10 # M
n <- 8 # k
Ex <- n*M/N
Ex
## [1] 1.6

Hypergeometric distribution - gene mutation (s13)

In a mutation test of 12 genes, 8 genes were found mutated in a cancer cell. If one selects 6 genes randomly, what is the expected number of mutated genes and the probability of finding 4 mutated genes out of the 6 genes? What if the 120 genes were tested initially and got 80 mutated genes?

N1 <- 12 # n = N-M
M1 <- 8 # M
n1 <- 6 # k
x1 <- 4 # x
Ex1 <- n1*M1/N1; Ex1
## [1] 4
# p(X=4)
dhyper(x1, M1, N1 - M1, n1)
## [1] 0.4545455
N2 <- 120 # n = N-M
M2 <- 80 # M
n2 <- 6 # k
x2 <- 4 # x

Ex2 <- n2*M2/N2; Ex2
## [1] 4
# p(X=4)
dhyper(x2, M2, N2 - M2, n2)
## [1] 0.3377274
# cf) for 1200
N3 <- 1200 # n = N-M
M3 <- 800 # M
n3 <- 6 # k
x3 <- 4 # x

Ex3 <- n3*M3/N3; Ex3
## [1] 4
# p(X=4)
dhyper(x3, M3, N3 - M3, n3)
## [1] 0.3300438
par(mfrow = c(1,3))
plot(0:12, dhyper(0:12, M1, N1 - M1, n1), type = 'h', 
     col = 3, ylim = c(0,0.5), xlab = 'x', ylab = 'f(x)',
     main = "Hypergeometric distribution (N=12)")
points(0:12, dhyper(0:12, M1, N1 - M1, n1), col = 2)
abline(h = 0, col = 3)

plot(0:12, dhyper(0:12, M2, N2 - M2, n2), type = 'h', 
     col = 3, ylim = c(0,0.5), xlab = 'x', ylab = 'f(x)',
     main = "Hypergeometric distribution (N=120)")
points(0:12, dhyper(0:12, M2, N2 - M2, n2), col = 2)
abline(h = 0, col = 3)

plot(0:12, dhyper(0:12, M3, N3 - M3, n3), type = 'h', 
     col = 3, ylim = c(0,0.5), xlab = 'x', ylab = 'f(x)',
     main = "Hypergeometric distribution (N=1200)")
points(0:12, dhyper(0:12, M3, N3 - M3, n3), col = 2)
abline(h = 0, col = 3)

par(mfrow = c(1,1))

dhyper(x1, M1, N1 - M1, n1); dhyper(x2, M2, N2 - M2, n2); dhyper(x3, M3, N3 - M3, n3)
## [1] 0.4545455
## [1] 0.3377274
## [1] 0.3300438

Negative Binomial distribution - baseball (s15)

The probability that Hanwha win against Doosan a game is 0.6, and Doosan win against Hanwha a game is 0.4. Find the probability that each team win the series by winning 4 games first out of 7 matches.

x <- 3 # number of failures x
size <- 4 # until y-th success denoted as size
prob <- 0.6 # probability of Hanwha win

# Hanhwa win
pnbinom(x, size, prob)
## [1] 0.710208
# Doosan win
pnbinom(x, size, 1 - prob)
## [1] 0.289792
par(mfrow = c(1,2))

plot(0:7, dnbinom(0:7, size, prob), type = 'h', 
     col = 3, ylim = c(0,0.3), xlab = 'x', ylab = 'f(x)',
     main = "Hanhwa win")
points(0:7, dnbinom(0:7, size, prob), col = 2)
abline(h = 0, col = 3)

plot(0:7, dnbinom(0:7, size, 1 - prob), type = 'h', 
     col = 3, ylim = c(0,0.3), xlab = 'x', ylab = 'f(x)',
     main = "Doosan win")
points(0:7, dnbinom(0:7, size, 1 - prob), col = 2)
abline(h = 0, col = 3)

par(mfrow = c(1,1))

Negative Binomial - gene expression (s15)

5% of the genes may express in a particular condition. What is the probability that we see two genes expressed under the same condition after examining at least five genes ? What is the mean trials to see two expressed genes?

x <- 0:2 # number of failures x
size <- 2  # until y-th success denoted as size
prob <- 0.05

# P(X >= 5) = 1-P(X<5) = 1 - {P(X=2) + P(X=3) + P(X=4)}
# P(X=0), P(X=1) is not considered due to the constraint X>2. 
1 - sum(dnbinom(x, size, prob))
## [1] 0.9859813
# mean trials
size/prob
## [1] 40

Geometric distribution - penalty kick (s17)

Assume that there is one soccer player, who made 20% of his attempted penalty kick goals in his career. Find the probability that he success first penalty kick after third attempts.

x <- 0:2 # number of failures (or number of trial - 1)
prob <- 0.2 # probability of success 

# P(X >= 3) = 1 - P(X < 4)
1 - sum(dgeom(x, prob))
## [1] 0.512

Geometric distribution - Cell division (s17)

Assume that the probability of producing mutation at each cell division is 0.10. What is the probability to produce a mutation at the third division, and average number of required cell division before getting the first mutation?

x <- 2 # number of failures (or number of trial - 1)
prob <- 0.1 # probability of success 

# P(X=3)
dgeom(x, prob)
## [1] 0.081
# E(x)
Ex <- 1/prob
Ex
## [1] 10
# Plot
plot(0:20, dgeom(0:20, prob), type = 'h', col = 3, ylim = c(0,0.15), 
     xlab = 'x', ylab = 'f(x)', main = "Geometric Distribution (p=0.1)")
points(0:20, dgeom(0:20, prob), col = 2)
abline(h = 0, col = 3)

Poisson distribution - software error (s21)

Suppose that the number of errors in a piece of software has a Poisson distribution with parameter \(\lambda = 3\). Find the probability that a piece of software has no errors and there are three or more errors in a piece of software.

# no error
dpois(0, lambda = 3)
## [1] 0.04978707
# three or more errors = 1 - less than 3 erros
1 - ppois(2, lambda = 3)
## [1] 0.5768099

Poisson distribution - gene mutation (s21)

Suppose that the number of mutation for 10 genes is expected to be 18 in a particular species. What is the probability of finding no mutation in a gene of this species?

# no mutation, lambda = 18 mutations / 10 genes
dpois(0, lambda = 1.8)
## [1] 0.1652989
# Plot
plot(0:10, dpois(0:10, lambda = 1.8), type = 'h', 
     col = 3, ylim = c(0,0.4), xlab = 'x', ylab = 'f(x)',
     main = "Poisson Distribution (p=1.8)")
points(0:10, dpois(0:10, lambda = 1.8), col = 2)
abline(h = 0, col = 3)

# Think about how the problem of Binomial Distribution can be transformed to
# that of Poisson Distribution.

n <- c(10, 100, 1000)
p <- c(0.18, 0.018, 0.0018)

dbinom(0, n, p)
## [1] 0.1374480 0.1626106 0.1650310
# approximating binomial distribution using poisson distribution 
# (when n is large, p is small)
lambda <- n*p
dpois(0, lambda)
## [1] 0.1652989 0.1652989 0.1652989

Exponential Distribution (s23)

Cell division occurs according to the Poisson process with a mean rate of two divisions per six minutes. What is the probability that it will take at least five minutes for the first division of the cell to take place?

  • \(X\): the waiting time in minutes until the first cell division take place

  • \(X\) has an exponential distribution with \(\lambda=2/6=1/3\)

# P(X > 5)
x <- 5
lambda <- 1/3

# f(x) = lambda*e^{-lambda * x}
1 - pexp(x, rate = lambda)
## [1] 0.1888756
x <- seq(0, 10, length = 100)
hx <- dexp(x, rate = lambda)
plot(x, hx, type = "l", lty = 2, xlab = "x value", ylim = c(0,0.5),
     ylab = "Density", main = "Exponential distribution (lambda = 1/3)")

Geometric Approximation of Exponential Distribution (s23)

Assume that the probability of producing mutation at each cell division in the question of slide 17 (Geometric Distribution) follows the Poisson process. Recalculate it by applying Exponential distribution. Assume that the cell division time is 10 minutes. Calculate the probability to produce a mutation within 30 min, and average waiting time in minutes until the first mutation?

# Exponential distribution
lambda <- 1/10
dexp(x = 2, rate = lambda)
## [1] 0.08187308
Ex <- 1/lambda
Ex
## [1] 10
# recall s17, Geometric approximation of exponential distribution
dgeom(x = 2, prob = 0.1)
## [1] 0.081
Ex <- 1/0.1
Ex
## [1] 10
# Plot
plot(0:20, dexp(0:20, lambda), type = 'h', col = 3, 
     ylim = c(0,0.15), xlab = 'x', ylab = 'f(x)', 
     main = "Approximation of exponential distribution by geometric distribution")
points(0:20, dexp(0:20, lambda), col = 2)
abline(h = 0, col = 3)
lines(0:20, dgeom(0:20, prob = 0.1), col = 1, lty = 3)
legend("topright", legend = c("Exponential", "Geometric"), col = c("black", "green"),
       lty = c(3, 1))

Gamma Distribution (s25)

Suppose we are interested in how long it will take for 20 metal sheets to be delivered to the panel construction lines where the delivery rate is 1.6 sheet/unit time. Find the expectation and variance of the waiting time.

\(f(x)= \frac{1}{(s^{a} \Gamma(a))} x^{(a-1)} e^{-(x/s)}\)

k <- 20 # shape = a = k
lambda <- 1.6 # scale = s = 1/lambda

Ex <- k/lambda
Var <- k/lambda^2
Ex; Var
## [1] 12.5
## [1] 7.8125
x <- seq(0, 25, length = 100)
hx <- dgamma(x, shape = k, scale = 1/lambda)
plot(x, hx, type = "l", lty = 2, xlab = "x value", ylim = c(0,0.2),
     ylab = "Density", main = "gamma distribution (a = 20, b = 0.625)")

Gamma approximation (s25)

Think about how the problem of Negative Binomial Distribution can be transformed to that of Gamma Distribution.

  • The Poisson process can be derived from the Binomial process by making n extremely large while p becomes very small, but within the constraint that np remains finite. In a Poisson process, the Gamma(a,b) distribution models the ‘time’ until observing a events where b is the mean time between events. The Negative binomial distribution is the binomial equivalent, modeling the number of failures to achieve s successes where \([(\frac 1 p)-1]\) is the mean number of failures per success. The Negative Binomial excludes the \(s\) successes which in terms of a Poisson process are not included in the waiting time because each event is assumed to be instantaneous. To make the two approaches exactly comparable, we should therefore think of the mean number of trials per success, equal to \((\frac 1 p)\). Then, we can make the following approximation:
# Gamma distribution
s <- 1
p <- 0.1 # lambda = 10
x <- seq(0, 100, length = 100)
hx <- dgamma(x, shape = s, scale = 1/p)

# Approximation of gamma by NB
plot(0:100, dnbinom(0:100, s, p), type = 'h', col = 3, 
     ylim = c(0,0.11), xlab = 'x', ylab = 'f(x)',
     main = "Approximation of Gamma distribution by Negative Binomial distribution")
points(0:100, dnbinom(0:100, s, p), col = 2)
abline(h = 0, col = 3)
lines(x, hx, lty = 3, col = 1)
legend("topright", legend = c("Gamma", "Negative Binomial"), 
       col = c("black", "green"), lty = c(3, 1))

for mathematical derivation, see this page

Beta Distribution - bee (s27)

When a queen bee leaves a bee colony to start a new hive, a certain proportion of the worker bees take flight and follow her. The proportion X of the worker bees that leave with the queen can be modeled by using a beta distribution with parameters a = 2.0 and b = 4.8. Then what could be the expected proportion and variance of bees leaving?

a <- 2.0
b <- 4.8
Ex <- a/(a + b)
Var <- a*b/((a + b)^2 * (a + b + 1))
Ex; Var
## [1] 0.2941176
## [1] 0.02661698
x <- seq(0, 1, length = 100)
hx <- dbeta(x, shape1 = a, shape2 = b)
plot(x, hx, type = "l", lty = 2, xlab = "x", ylim = c(0,2.5),
     ylab = "f(x)", main = "beta distribution (a = 2.0, b = 4.8)")

# What is the probability that more than 29% of worker bee leave with queen?
pbeta(0.29, shape1 = a, shape2 = b)
## [1] 0.5378305

Bayesian Estimation (s27)

Let’s solve the above problem in Bayesian way. Beta distributions can be used to model the prior probability distribution of events. The beta distribution is the conjugate prior of the binomial distribution, and our problem is the binomial distribution problem, so the posterior distribution is also the beta distribution. In fact, the \(\alpha\) and \(\beta\) parameters of the beta distribution can be viewed as pseudocounts, which are added to the observed counts to reflect prior knowledge.(for more information about conjugate prior, see this page)

# beta prior
alpha <- 2.0
beta <- 4.8
plot(seq(0,1,by = 0.01), dbeta(seq(0,1,by = 0.01), alpha, beta), type = "l",
     ylim = c(0, 5), xlim = c(0,1), ylab = "Density", xlab = "x")

# set the observed number of successes and failures.
n.success = 3  # number of bees follow the queen
n.failure = 7  # number of bees did not follow the queen

a = alpha + n.success  # convert to beta parameter
b = beta + n.failure   # convert to beta parameter

# with weak observation
theta <- seq(0.005, 0.995, length = 500)
plot(theta, dbeta(theta, a, b), type = "l", col = "red",
     ylim = c(0, 5), xlim = c(0,1), ylab = "Density", xlab = "x")
lines(theta, dbeta(theta, n.success + 1, n.failure + 1), 
      lty = 2, col = "green")
lines(theta, dbeta(theta, alpha, beta), lty = 3, col = "blue")
legend("topright", legend = c("Prior", "likelihood", "Posterior"), 
       col = c("blue", "green", "red"), lty = c(3, 2, 1))

# Add more observations
n.success = 30  # number of bees follow the queen
n.failure = 70  # number of bees did not follow the queen

a = alpha + n.success  # update parameter
b = beta + n.failure   # update parameter

# with strong observation
plot(theta, dbeta(theta, a, b), type = "l", col = "red",
     ylim = c(0, 10), xlim = c(0,1), ylab = "Density", xlab = "x")
lines(theta, dbeta(theta, n.success + 1, n.failure + 1), 
      lty = 2, col = "green")
lines(theta, dbeta(theta, alpha, beta), lty = 3, col = "blue")
legend("topleft", legend = c("Prior", "likelihood", "Posterior"), 
       col = c("blue", "green", "red"), lty = c(3, 2, 1))

In conclusion, The Bayes estimate in this example is simply the mean of beta distribution \(E[\theta]\), which should be \(\theta = \frac{a+X}{b+n-X}\). for more information about this topic, plz refer to this link

Bayesian Estimation using package - Advanced learner only (s27)

You can also use functions which provided by LearnBayes package

# install.packages("LearnBayes")
library(LearnBayes)

## calculate some properties of the Beta distribution
calcBetaMode = function(aa, bb) {
  beta.mode = (aa - 1)/(aa + bb - 2)
  return(beta.mode)
}
calcBetaMean = function(aa, bb) {
  beta.mean = (aa)/(aa + bb)
  return(beta.mean)
}
calcBetaVar = function(aa, bb) {
  beta.var = (aa * bb)/(((aa + bb)^2) * (aa + bb + 1))
  return(beta.var)
}
calcBetaMedian = function(aa, bb) {
  beta.med = (aa - 1/3)/(aa + bb - 2/3)
  return(beta.med)
}
calcBetaSkew = function(aa, bb) {
  beta.skew = (2 * (bb - aa) * sqrt(aa + bb + 1) ) / ((aa + bb + 2)/sqrt(aa + bb))
  return(beta.skew)
}


priorToPosterior = function(successes, total, a, b) {
  ## Note the rule of succession
  likelihood.a = successes + 1
  likelihood.b = total - successes + 1
  
  ## Create posterior
  posterior.a = a + successes;
  posterior.b = b + total - successes
  theta = seq(0.005, 0.995, length = 500)
  
  ## Calc density
  prior = dbeta(theta, a, b)
  likelihood = dbeta(theta, likelihood.a, likelihood.b)
  posterior = dbeta(theta, posterior.a, posterior.b)
  
  ## Plot prior, likelihood, and posterior
  
  ## Can be used to scale down the graph if desired.
  ## However, the density is different for each prior, likelihood, posterior
  m.orig = apply( cbind(prior, likelihood, posterior), 2, max)
  m = max(c(prior, likelihood, posterior))
  
  plot(theta, posterior, type = "l", ylab = "Density", lty = 2, lwd = 3,
       main = paste("Prior: beta(", round(a,2), ",", round(b,2), "); Data: B(", total, ",", successes, "); ",
                    "Posterior: beta(", round(posterior.a,2), ",", round(posterior.b,2), ")", sep = ""), ylim = c(0, m), col = 1)
  lines(theta, likelihood, lty = 1, lwd = 3, col = 2)
  lines(theta, prior, lty = 3, lwd = 3, col = 3)
  legend("topleft",y = m, c("Prior", "Likelihood", "Posterior"), lty = c(3, 1, 2),
         lwd = c(3, 3, 3), col = c(3, 2, 1))
  
  prior.mode = calcBetaMode(a, b)
  likelihood.mode = calcBetaMode(likelihood.a, likelihood.b)
  posterior.mode = calcBetaMode(posterior.a, posterior.b)
  prior.mean = calcBetaMean(a, b)
  likelihood.mean = calcBetaMean(likelihood.a, likelihood.b)
  posterior.mean = calcBetaMean(posterior.a, posterior.b)
  prior.med = calcBetaMedian(a, b)
  likelihood.med = calcBetaMedian(likelihood.a, likelihood.b)
  posterior.med = calcBetaMedian(posterior.a, posterior.b)
  prior.var = calcBetaVar(a, b)
  likelihood.var = calcBetaVar(likelihood.a, likelihood.b)
  posterior.var = calcBetaVar(posterior.a, posterior.b)
  prior.skew = calcBetaSkew(a, b)
  likelihood.skew = calcBetaSkew(likelihood.a, likelihood.b)
  posterior.skew = calcBetaSkew(posterior.a, posterior.b)
  
  print(paste("Mode: prior=",prior.mode,"; Likelihood=",likelihood.mode,"; Posterior=",posterior.mode))
  print(paste("Mean: prior=",prior.mean,"; Likelihood=",likelihood.mean,"; Posterior=",posterior.mean))
  print(paste("~Approx Median (for a and b > 1): prior=",prior.med,"; Likelihood=",likelihood.med,", for Posterior=",posterior.med))
  print(paste("Var: prior=",prior.var,"; Likelihood=", likelihood.var,"; Posterior=",posterior.var))
  print(paste("Skewness: prior=",prior.skew,"; Likelihood=",likelihood.skew,"; Posterior=",posterior.skew))
  return(list(a = posterior.a,b = posterior.b))
}

# calculate posterior for proportion
calcPosteriorForProportion <- function(successes, total, a, b)
{
  # Adapted from triplot() in the LearnBayes package
  # Plot the prior, likelihood and posterior:
  likelihood_a = successes + 1; likelihood_b = total - successes + 1
  posterior_a = a + successes;  posterior_b = b + total - successes
  theta = seq(0.005, 0.995, length = 500)
  prior = dbeta(theta, a, b)
  likelihood = dbeta(theta, likelihood_a, likelihood_b)
  posterior  = dbeta(theta, posterior_a, posterior_b)
  m = max(c(prior, likelihood, posterior))
  plot(theta, posterior, type = "l", ylab = "Density", lty = 2, lwd = 3,
       main = paste("beta(", a, ",", b, ") prior, B(", total, ",", successes, ") data,",
                    "beta(", posterior_a, ",", posterior_b, ") posterior"), ylim = c(0, m), col = "red")
  lines(theta, likelihood, lty = 1, lwd = 3, col = "blue")
  lines(theta, prior, lty = 3, lwd = 3, col = "green")
  legend(x = 0.8, y = m, c("Prior", "Likelihood", "Posterior"), lty = c(3, 1, 2),
         lwd = c(3, 3, 3), col = c("green", "blue", "red"))
  # Print out summary statistics for the prior, likelihood and posterior:
  calcBetaMode <- function(aa, bb) { BetaMode <- (aa - 1)/(aa + bb - 2); return(BetaMode); }
  calcBetaMean <- function(aa, bb) { BetaMean <- (aa)/(aa + bb); return(BetaMean); }
  calcBetaSd   <- function(aa, bb) { BetaSd <- sqrt((aa * bb)/(((aa + bb)^2) * (aa + bb + 1))); return(BetaSd); }
  prior_mode      <- calcBetaMode(a, b)
  likelihood_mode <- calcBetaMode(likelihood_a, likelihood_b)
  posterior_mode  <- calcBetaMode(posterior_a, posterior_b)
  prior_mean      <- calcBetaMean(a, b)
  likelihood_mean <- calcBetaMean(likelihood_a, likelihood_b)
  posterior_mean  <- calcBetaMean(posterior_a, posterior_b)
  prior_sd        <- calcBetaSd(a, b)
  likelihood_sd   <- calcBetaSd(likelihood_a, likelihood_b)
  posterior_sd    <- calcBetaSd(posterior_a, posterior_b)
  print(paste("mode for prior=",prior_mode,", for likelihood=",likelihood_mode,", for posterior=",posterior_mode))
  print(paste("mean for prior=",prior_mean,", for likelihood=",likelihood_mean,", for posterior=",posterior_mean))
  print(paste("sd for prior=",prior_sd,", for likelihood=",likelihood_sd,", for posterior=",posterior_sd))
}

# weak observation
n.success = 3  # number of bees follow the queen
n.failure = 7  # number of bees did not follow the queen

posterior.out = priorToPosterior(n.success, n.success + n.failure, alpha, beta)

## [1] "Mode: prior= 0.208333333333333 ; Likelihood= 0.3 ; Posterior= 0.27027027027027"
## [1] "Mean: prior= 0.294117647058824 ; Likelihood= 0.333333333333333 ; Posterior= 0.297619047619048"
## [1] "~Approx Median (for a and b > 1): prior= 0.271739130434783 ; Likelihood= 0.323529411764706 , for Posterior= 0.289256198347107"
## [1] "Var: prior= 0.0266169816342827 ; Likelihood= 0.0170940170940171 ; Posterior= 0.0117439297816505"
## [1] "Skewness: prior= 4.63454509789934 ; Likelihood= 7.13714056959817 ; Posterior= 12.5096656156832"
# Add more observations
n.success = 30  # number of bees follow the queen
n.failure = 70  # number of bees did not follow the queen

posterior.out = priorToPosterior(n.success, n.success + n.failure, alpha, beta)

## [1] "Mode: prior= 0.208333333333333 ; Likelihood= 0.3 ; Posterior= 0.295801526717557"
## [1] "Mean: prior= 0.294117647058824 ; Likelihood= 0.303921568627451 ; Posterior= 0.299625468164794"
## [1] "~Approx Median (for a and b > 1): prior= 0.271739130434783 ; Likelihood= 0.302631578947368 , for Posterior= 0.298366834170854"
## [1] "Var: prior= 0.0266169816342827 ; Likelihood= 0.00205391503641243 ; Posterior= 0.00194666091829148"
## [1] "Skewness: prior= 4.63454509789934 ; Likelihood= 78.8452157542693 ; Posterior= 84.4189363883143"