Jim — Sep 22, 2013, 5:27 PM
# MATH 6480 - Fall 2013
# Learning about a Normal Distribution
# Suppose sleeping times of college students are N(mu, sigma).
# We assume that priors on mu and sigma^2 are independent where
# mu is distributed N(mu0, sigma0) and sigma^2 is distributed IG(a, b)
# Here are some prior beliefs.
# I am 90% confident that mu is between 7 and 9 hours.
# I am 90% confident that sigma is between 1.5 and 3 hours.
# First, use LearnBayes function norm.select to find matching normal prior
library(LearnBayes)
normal.select(list(p=.05, x=7), list(p=.95, x=9))
$mu
[1] 8
$sigma
[1] 0.608
# Matching normal prior is N(8, 0.61)
# Next, want to find a matching Inverse Gamma prior that matches beliefs.
# I'll use a roundabout way of figuring out a IG prior which approximately
# matches these beliefs.
# Suppose sigma is N(mu.s, sigma.s)
normal.select(list(p=.05, x=1.5), list(p=.95, x=3.5))
$mu
[1] 2.5
$sigma
[1] 0.608
# See that mu.s = 2.5, sigma.s = 0.61
# Simulate 1000 values from this normal distribution.
# Square values to get 1000 values from prior on precision P = 1/sigma^2.
# Fit this "data" to a gamma distribution to find matching a, b.
P <- 1 / rnorm(10000, 2.5, 0.61) ^ 2
library(MASS)
fitdistr(P, dgamma, list(shape=1, rate=1))
Warning: NaNs produced Warning: NaNs produced Warning: NaNs produced
Warning: NaNs produced Warning: NaNs produced Warning: NaNs produced
shape rate
2.39240 11.10183
( 0.03176) ( 0.16392)
# Matching prior on sigma^2 is IG(2.8, 13.5)
# 5th and 9th quantiles of this prior on sigma are 1.5 and 4.3
# This will be fine for our purposes.
# Now we observe sleeping data from the studentdata data frame in
# LearnBayes.
y <- with(studentdata, WakeUp - ToSleep)
# We remove some missing values from the data.
y <- y[complete.cases(y)]
# To check if it is normally distributed, construct a normal
# probability plot.
qqnorm(y); qqline(y)
# Although we notice the discrete nature of the data, normality
# seems to be a reasonable assumption.
# Simulate from posterior distribution of (mu, sigma2) using
# Gibbs sampling.
# The function normpostsim will simulate values from the posterior
# of (mu, sigma2) using Gibbs sampling.
S <- normpostsim(y, prior=list(mu=c(8.0, 0.61), sigma2=c(2.8, 13.5)))
# Construct a scatterplot of the simulated draws of (mu, sigma^2).
plot(S$mu, S$sigma2)
