Instruction: copy and paste the following codes in R first. Or save the whole code
in a file, e.g. poisson.r, then source it in R: > source('poisson.r')

    R-codes for Poisson sample with conjugate priors for the mean rate
    R-codes for Binomial sample with conjugate priors for the success proportion
    R-codes for Normal sample with conjugate priors for the mean
    R-codes for Exponential sample with conjugate priors for the mean rate

--------------------------------------------------------------------------------------

1. Class Example: Transmission Error Data

Data: 6 one-hour observation periods with 1, 0, 1, 2, 1, 0 transmission errors. The data can
be modeled as a sample from a Poisson distribution with mean rate lambda.
The data on previous system established the mean error rate as 1.6 errors per hour and the new
system has a design goal of cutting this rate in half to .8 errors per hour.
An expert has a prior belief that the median of the new system should be close to .8. Consider
a gamma distribution with alpha=2 and beta=1/2.1 as the prior distribution. Does the new system
seem to be able to achieve the design goal?

Likelihood: Poisson (lambda)
> tran.dat <- c(1, 0, 1, 2, 1, 0)

Conjugate prior: Gamma (2, 1/2.1)
Obtain the Posterior distribution with the mean, variance, and the mode:

# To obtain the conjugate posterior distribution, you need to specify only data set and parameters
in the function 'pg.sum'

>pg.sum( tran.dat,2,2.1)   #Note that second parameter in Gamma distn is specified as 1/beta in R.

# $alfa1 [1] 7 $beta1 [1] 8.1         # So the Posterior distribution is Gamma(7, 1/8.1)
# $Mean [1] 0.8641975 $Var [1] 0.1066911 $Moodi [1] 0.7407407      # Posterior inference

Triplot of the prior, likelihood, and the posterior distribution in one figure:
> pg.trip( tran.dat, 2, 2.1)

Or you can plot them on your own:

> lam <-seq (0, 3, length=50 )
> prior.g <-dgamma(lam ,2, 2.1)
> post.g<-dgamma(lam, 7, 8.1)

> plot(lam, post.g, type="l",col="1" )
> lines(lam, prior.g, type="l", col="2")

Posterior probability :
>  pgamma(0.8,  7 , 8.1)        # P( meet the design goal)
> 1- pgamma(1.6,  7 , 8.1)     # P( worse than the old system)

90% Posterior highest density (credible) Interval:
> pg.hdi( tran.dat, 2, 2.1, .9)

-----------------------------------------------------------------------------------------
2. Example: Clinical trial data
A new treatment protocol is proposed for a particular form of cancer. The measure of
success will be the proportion of patients that survive longer than six months after
diagnosis. With the present treatment, this success rate is 40%. Letting theta be the
success rate of the new treatment, a doctor assesses her prior beliefs about theta as
follows. She judges that her expectation of theta is E(theta)=0.45, and her standard
deviation is 0.07. Assume that her beliefs can be represented by a beta distribution.

A clinical trial of the new treatment protocol is carried out. Out of 70 patients
in the trial, 34 survive beyond six months from diagnosis.Should the new treatment
protocol replace the old one?

Likelihood:  Binomial (n, theta)
> clin.dat <- c(rep(1, 34), rep(0,36) )

Conjugate prior for theta:  Beta(22.28, 27.23)
Obtain the Posterior distribution with the mean, variance, and the mode:

> bb.sum(clin.dat, 1, 22.28,27.23)
#$alfa1[1] 56.28 $beta1[1] 63.23    #  i.e. Posterior: Beta (56.28, 63.23)
#$Mean[1] 0.4709229  $Var[1] 0.002067501 $Moodi [1] 0.4704280

-----------------------------------------------------------------------------------------

3. Example:  Earth's density ( Henry Cavendish )

Likelihood: N( mu, 0.25)  # Data

> den.dat <- c( 5.36, 5.29, 5.58, 5.65, 5.57, 5.53, 5.62, 5.29, 5.44, 5.34, 5.79,
            5.10, 5.27,  5.39, 5.42, 5.47, 5.63, 5.34, 5.46, 5.30, 5.78, 5.68, 5.85)
> summary(den.dat)
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  5.100   5.340   5.460   5.485   5.625   5.850

Conjugate prior for mu:  N(5,4)

Obtain the Posterior distribution with the mean, variance, and the mode:
> nn.sum(den.dat, 5,4, 0.25)
#$Mean [1] 5.483469 $Var[1]0.010840  i.e. Posterior: Normal(5.4835, 0.01084)
 

Plot of the prior, likelihood, and the posterior distribution

> x <- seq(0,10,50)
> x.prior <- dnorm(x, 5, 2)    # 2= std.dev :  sqrt(4)
> x.post <- dnorm(x, 5.486, sqrt(0.01084) )
> x.lkhd <-  dnorm(x, mean(x),sqrt(0.25/ length(x)) ) # mean=x_bar, variance=sigma/n

> plot(x, x.prior, ylim=c(0,4), type='l', col=4)
> lines(x, x.post, lty=2, col=2)
> lines(x, x.lkhd,  lty=3, col=3)
> legend(0, 3, c('prior','posterior','lkhd'), col=c(4,2,3),  lty=1:3 )

Triplot of the prior, likelihood, and the posterior distribution in one figure using the function 'nn.trip':
> nn.trip(eden.dat, 5,4, 0.25)

90% Posterior highest density (credible) Interval:
> nn.hdi(eden.dat, 5,4, 0.25, .90)

Posterior  probability:
> 1-pnorm(5.5, 5.4835, sqrt(0.01084) )
#[1] 0.43704
________________________________________________________________________________

Your turn

1. A teacher administers a standardized aptitude exam to a class of 20 students.
She knows that the national distribution of scores for this exam is normal with mean 65
and standard deviation 9. She is willing to assume that students in her class have
aptitude values that are normally distributed with standard deviation 9, but with
unknown mean.

Her students come from an affluent, educated population. She assesses a 95% credible
(highest density) interval for the mean of her class distribution as (60,80). She
therefore assumes a priori that the mean of her class distribution is normal with mean 70
and standard deviation 5.

When the scores are returned, she opens the first 5 envelopes and sees scores of
78,80,95,81, and 77.

i) Find the posterior distribution of class mean.

ii) Obtain 95% highest posterior density interval for the mean.

iii) Make three plots of prior, likelihood, and posterior distributions in one figure.
Give some comments on these plots.

iv) Give a predictive distribution for the next observation using the code 'nn.pred'.
 

2. The following data are failure times measured in CPU seconds for a real-time command
and control system. That is, the first failure occurred 3 seconds after the start of the
test, the second 30 seconds later at 33 seconds, ... the last occurred at 759 seconds.
The data consists of the following number of seconds between errors:

   3, 30, 113, 81, 115, 9, 2, 91, 112, 15, 138, 50

Assume that the times between failures are independent and distributed exponentially
with a constant failure rate. Assume a Gamma(1,1) prior distribution.

i) Note that the Gamma distribution is the conjugate prior for the exponential
distribution of the unknown failure rate. First, obtain the analytic form of the
posterior distribution using conjugacy by hand.

 p(x_i|theta) = theta* exp(-x_i * theta),      theta ~exp(theta)

ii) Obtain the posterior distribution using the code given above and fill in the blanks below.

                               Mean         Standard Deviation
   Prior
   Posterior

iii)  Set the data x as the sequence of 50 values from 0 to 0.10 first.
Obtain the probability density functions for the prior and posterior distributions.

iv) Make plots of prior and posterior distributions overlaid in one figure.
Comment on the plots, including the central tendency of each of the plotted functions and
the amount of spread in each of the plotted function.

v) The software test continued for 91208 CPU seconds. The below shows a sequence
of 13 additional observations taken later during the test. Use these observations to
create a data set of 12 times between failures.

26, 114, 325, 55, 242, 68, 422, 180, 10, 1146, 600, 15

Under the assumption of this problem that the times between failures are independent,
exponentially distributed with constant failure rate, it would be appropriate to use the
posterior distribution from i) (or ii)) as your prior distribution for this data set.

Obtain the posterior distribution and make two plots as in iv). Comment on your results.

vi) Using a Gamma(1,1) prior distribution, repeat the problem iv). Comment on your
results.

3. The numbers of misprints (typos) spotted on the first few pages of an early draft of a
book were

3, 4, 2, 1, 2, 3

It seems reasonable that these numbers should constitute a sample from a Poisson
distribution of unknown mean lambda. From the past experience of typing such a book,
the author would like to use a prior with a mean about 3 and variance about 2.

i) Using the conjugate prior distribution, find the posterior distribution.

ii) Give the plots of prior and the posterior distribution on one page.

4. An usual condition of pregnancy called placenta previa obstructs normal baby delivery.
An early study concerning the gender (sex) of placenta previa births in Germany found
that a total of 980 births, 437 were female.

i) Obtain the posterior distribution using a conjugate prior with alpha 0.5 and beta 1/12.

ii) Find the 95% highest posterior density interval for the unknown proportion of placenta
previa births.

iii) How much evidence does this provide for the claim that the proportion of female births in
the population of placenta previa births is lessthan 0.485 (the currently assumed proportion
of female births in the general population)? Compute the appropriate posterior probability.