## Tuesday 7 March 2017
data = read.table('notes_models_poisproc.dat')
eventtimes = data[,1]
y = length(eventtimes)
T = 24*7
thetaMLE = y / T
theta = seq(0,2*thetaMLE,length.out=1000)
posterior = dgamma(theta,y,rate=T)
plot(theta,posterior,'l')
plot(eventtimes,1:y,pch=20,xaxp=c(0,24*7,7))
waittimes = eventtimes[2:y] - eventtimes[1:y-1]
library(rethinking)
dens(waittimes)
wt = seq(0,1,length.out=100)
pt = dexp(wt,thetaMLE)
lines(wt,pt,col='green')
plot(eventtimes[2:y],waittimes,pch=20,
    xaxp=c(0,24*7,7))
plot(eventtimes[2:y] %% 24,waittimes,pch=20,
    xaxp=c(0,24,6))
hist(eventtimes %% 24,breaks = 24)
hist(abs(eventtimes %% 24 - 12))
yday = sum(abs(eventtimes %% 24 - 12)<6)
yday
ynight = sum(abs(eventtimes %% 24 - 12)>6)
ynight
sum(eventtimes %% 24 == 6)
sum(eventtimes %% 24 == 18)
yday = yday + sum(eventtimes %% 24 == 6)
ynight = ynight + sum(eventtimes %% 24 == 18)
yday
ynight
Tday = 12*7
Tnight = 12*7
postday = dgamma(theta,yday,rate=Tday)
postnight = dgamma(theta,ynight,rate=Tnight)
plot(theta,postnight,'l',col='blue')
lines(theta,postday,col='red')
pbinom(ynight,y,0.5)
logB1 = lgamma(yday) + lgamma(ynight) - lgamma(y)
logB1
lbeta(yday,ynight)
logB2 = ( y * log(T)
  - yday * log(Tday) - ynight * log(Tnight) )
logB2
y * log(2)
logOccam = -log(log(1e4))
logOccam
logB = logB1 + logB2 + logOccam
logB
exp(logB)
## Thursday 9 March 2017
library(rethinking)
data(Howell1)
y = Howell1[Howell1$age >= 18,'weight']
ybar = mean(y)
ybar
sy = sd(y)
sy
n = length(y)
mu = seq(ybar-5*sy/sqrt(n),ybar+5*sy/sqrt(n),
   length.out=1000)
t = (mu-ybar)/(sy/sqrt(n))
pmu = (sqrt(n)/sy) * dt(t,n-1)
pgauss = dnorm(mu,mean=ybar,sd=sy/sqrt(n))
plot(mu,pmu,'l',col='black')
lines(mu,pgauss,col='blue')
hist(y,breaks=sqrt(length(y)))
x = Howell1[Howell1$age >= 18,'height']
plot(x,y,'p',pch=20)
xbar = mean(x)
Sxx = sum((x-xbar)^2)
Sxy = sum((x-xbar)*(y-ybar))
Syy = sum((y-ybar)^2)
s = sqrt((Syy-Sxy^2/Sxx)/(n-2))
print(s)
print(sy)
alpha = mu
palpha = dnorm(alpha,mean=ybar,sd=s/sqrt(n))
plot(alpha,palpha,'l',col='brown')
lines(mu,pmu,col='blue')
betahat = Sxy/Sxx
print(betahat)
sigbeta = s / sqrt(Sxx)
print(sigbeta)
betavals = seq(0,betahat+5*sigbeta,length.out=1000)
pbeta = dnorm(betavals,mean=betahat,sd=sigbeta)
plot(betavals,pbeta,'l')