## Thursday 26 January 2017
theta = seq(from=0, to=1, length.out = 1000)
d_theta = theta[2]-theta[1]
prior = dunif(theta)
ytot = 3
n = 4
likelihood = dbinom(ytot,n,theta)
posterior = prior * likelihood
posterior = posterior / (sum(posterior) * d_theta)
plot(theta,posterior,'l')
sum(theta*posterior)*d_theta
## Tuesday 7 February 2017
n = 4
ytot = 3
N = 5
thetasample = rbeta(N, ytot+1, n-ytot+1)
thetasample
thetasample = rbeta(N, ytot+1, n-ytot+1)
thetasample
set.seed(20170207)
thetasample = rbeta(N, ytot+1, n-ytot+1)
thetasample
set.seed(20170207)
thetasample = rbeta(N, ytot+1, n-ytot+1)
thetasample
set.seed(20170206)
thetasample = rbeta(N, ytot+1, n-ytot+1)
thetasample
set.seed(20170207)
N = 10000
thetasample = rbeta(N, ytot+1, n-ytot+1)
plot(thetasample,type='p',pch='.')
hist(thetasample)
library(rethinking)
dens(thetasample)
d_theta = 0.001
theta = seq(from=0, to=1, length.out=1000)
posterior = dbeta(theta,ytot+1,n-ytot+1)
lines(theta,posterior,col="blue")
posterior2 = dbinom(ytot,n,theta)
posterior2 = posterior2/sum(posterior2)
plot(theta,posterior2,type='l')
thetasample2 = sample(theta,size=N,prob=posterior2,
                      replace=TRUE)
dens(thetasample2)
lines(theta,posterior2/d_theta,col="blue")
mean(thetasample2)
sum(posterior2*theta)
median(thetasample2)
var(thetasample2)
sum(posterior2*(theta-sum(posterior2*theta))^2)
sd(thetasample2)
mean(thetasample2<0.5)
sum(posterior2[theta<0.5])
mean(thetasample);mean(thetasample2)
median(thetasample);median(thetasample2)
sd(thetasample);sd(thetasample2)
mean(thetasample<0.5)
HPDI(thetasample2,prob=0.95)
dens(thetasample2,show.HPDI = 0.95)
lambdasample=logit(thetasample)
dens(lambdasample)
lambda3 = seq(from=-5,to=5,length.out=1000)
pdf_lambda3 = (
   (1+exp(-lambda3))^(-(ytot+1))
   * (1+exp(lambda3))^(-(n-ytot+1))
 ) / ( beta(ytot+1,n-ytot+1) )
lines(lambda3,pdf_lambda3,col="blue")
## Thursday 9 February 2017
H = matrix(c(3,-2,-2,2),ncol=2)
theta1hat = 4
theta2hat = 3
Ngrid = 101
theta1 = seq(from=theta1hat-4,to=theta1hat+4,
    length.out=Ngrid)
theta2 = seq(from=theta2hat-4,to=theta2hat+4,
    length.out=Ngrid)
ones = rep(1,Ngrid)
theta1grid = outer(theta1,ones,'*')
theta2grid = outer(ones,theta2,'*')
theta1grid[1:5,1:5]
theta2grid[1:5,1:5]
logpost = -0.5 * (
    H[1,1] * (theta1grid-theta1hat)^2
    + 2 * H[1,2] * (theta1grid-theta1hat)
        *(theta2grid-theta2hat)
    + H[2,2] * (theta2grid-theta2hat)^2
)
max(logpost)
logpost = logpost - max(logpost)
contour(theta1,theta2,logpost,levels=0:-10)
posterior = exp(logpost)
posterior = posterior / sum(posterior)
contour(theta1,theta2,posterior)
post_marg1 = 0.0 * theta1
for (i in 1:Ngrid){post_marg1[i]=sum(posterior[i,])}
post_max1 = posterior[theta2grid==theta2hat]
post_max1 = post_max1/sum(post_max1)
plot(theta1,post_max1,type='l',col='red')
lines(theta1,post_marg1,col='blue')
sum(posterior*(theta1grid-theta1hat)^2)
sum(post_marg1*(theta1-theta1hat)^2)
sum(post_max1*(theta1-theta1hat)^2)
post_marg2 = 0.0 * theta2
for (i in 1:Ngrid){post_marg2[i]=sum(posterior[,i])}
post_max2 = posterior[theta1grid==theta1hat]
post_max2 = post_max2/sum(post_max2)
plot(theta2,post_max2,type='l',col='red')
lines(theta2,post_marg2,col='blue')
sum(posterior*(theta2grid-theta2hat)^2)
sum(post_marg2*(theta2-theta2hat)^2)
sum(post_max2*(theta2-theta2hat)^2)
sum(posterior*(theta1grid-theta1hat)
    *(theta2grid-theta2hat))
solve(H)
N = 10000
set.seed(20170209)
idxsample = sample.int(n=Ngrid*Ngrid,size=N,
    prob=posterior,replace=TRUE)
theta1sample = theta1grid[idxsample]
theta2sample = theta2grid[idxsample]
plot(theta1sample,theta2sample,pch='.')
dtheta1 = theta1[2]-theta1[1]
dtheta2 = theta2[2]-theta2[1]
theta1sample = ( theta1grid[idxsample]
    + runif(N,min=-0.5*dtheta1,
            max=0.5*dtheta1) )
theta2sample = ( theta2grid[idxsample]
    + runif(N,min=-0.5*dtheta2,
            max=0.5*dtheta2) )
plot(theta1sample,theta2sample,pch='.')
library(rethinking)
dens(theta1sample)
lines(theta1,post_max1/dtheta1,col='red')
lines(theta1,post_marg1/dtheta1,col='blue')
dens(theta2sample)
lines(theta2,post_max2/dtheta2,col='red')
lines(theta2,post_marg2/dtheta2,col='blue')
mean(theta1sample)
mean(theta2sample)
var(theta1sample)
var(theta2sample)
cov(theta1sample,theta2sample)
## Tuesday 14 February 2017
Ngrid = 101
t = seq(from=-6,to=6,length.out=Ngrid)
ones = rep(1,Ngrid)
tgrid = outer(t,ones,'*')
n = 5
nu = n-1
chisqmax = nu + 6*sqrt(2*nu)
chisqmin = nu*(nu/chisqmax)^1.4
logchisq = seq(from=log(chisqmin),to=log(chisqmax),
   length.out=Ngrid)
logchisqgrid = outer(ones,logchisq,'*')
chisqgrid = exp(logchisqgrid)
logpost = ( 0.5*n*logchisqgrid
	- 0.5 * (1+tgrid^2/nu) * chisqgrid
	)
max(logpost)
logpost = logpost - max(logpost)
posterior = exp(logpost)
posterior = posterior / sum(posterior)
contour(t,logchisq,posterior)
chisq = exp(logchisq)
contour(t,chisq,posterior/chisqgrid)
set.seed(20170214)
N=10000
idxsample = sample.int(n=Ngrid*Ngrid,size=N,
    prob=posterior,replace=TRUE)
dt = t[2] - t[1]
dlogchisq = logchisq[2] - logchisq[1]
tsample = ( tgrid[idxsample]
    + runif(N,min=-0.5*dt,
            max=0.5*dt) )
logchisqsample = ( logchisqgrid[idxsample]
    + runif(N,min=-0.5*dlogchisq,
            max=0.5*dlogchisq) )
plot(tsample,logchisqsample,pch='.')
chisqsample = exp(logchisqsample)
plot(tsample,chisqsample,pch='.')
## Thursday 16 February 2017
library(gtools)
alpha = c(1,1,1)
N = 2000
set.seed(20170216)
thetasample = rdirichlet(N,alpha)
library(ggtern)
myframe = function(mysample) {
   ( data.frame(theta1=mysample[,1],
        theta2=mysample[,2],
        theta3=mysample[,3]) )
}
tickvals = seq(.2,1,.2)
addmyopts = function(myplot) {
   ( myplot
     + geom_point(size=0.01)
     + tern_limits(breaks=tickvals,labels=tickvals)
     + xlab(expression(theta[1]))
     + ylab(expression(theta[2]))
     + zlab(expression(theta[3]))
     + theme_light()
)
}
thetaframe = myframe(rdirichlet(N,alpha))
addmyopts(ggtern(thetaframe,aes(theta1,theta2,theta3)))
alpha = c(0.1,0.1,0.1)
thetaframe = myframe(rdirichlet(N,alpha))
addmyopts(ggtern(thetaframe,aes(theta1,theta2,theta3)))
alpha = c(200,500,300)
thetaframe = myframe(rdirichlet(N,alpha))
addmyopts(ggtern(thetaframe,aes(theta1,theta2,theta3)))
alpha = c(20,50,30)
thetaframe = myframe(rdirichlet(N,alpha))
addmyopts(ggtern(thetaframe,aes(theta1,theta2,theta3)))
alpha = c(2,5,3)
thetaframe = myframe(rdirichlet(N,alpha))
addmyopts(ggtern(thetaframe,aes(theta1,theta2,theta3)))