# import the data
mar=read.delim("marmotbodymeasurements.txt")

# select a subset of the data: only adult animals
maradult=mar$Peso[mar$Age=="ADULTO"]

# the hist of the weight of the adults looks approx normal
hist(maradult,nclass=10,prob=T)

# compute the posterior mean
(32*mean(maradult)+20*3500)/52 = 3509.615

# plot the posterior distribution for sigma^(-2)
curve(dgamma(x,21,3645096),from=0,to=10^(-5))

# simulation from the marginal posterior of mu (it is a T-distribution)
sigma=1/rgamma(10^5,21,3645096)
mu=rnorm(10^5)*sqrt(sigma/52)+3509.615
hist(mu,nclass=200,col="wheat")

# if we change the strenght of the prior belief by going
# from 5 to 1  and from 20 to 2 (number of imaginary observations) we get ...
(sum((marmotadult-mean(marmotadult))^2)+1*(250)^2+32*2*(3515.625-3500)^2/(2+32))/2
sigmas=1/rgamma(10^5,17,3518824)
mus=rnorm(10^5)*sqrt(sigmas/34)+3514.706

# the predictive distribution of the weight of a new marmot
newmarmot=rnorm(10^5)*sqrt(sigmas)+mus
hist(newmarmot,nclass=200,col="chocolate")

# Comparison with older parameters
sigmas=1/rgamma(10^5,21,3645096)
mus=rnorm(10^5)*sqrt(sigmas/52)+3509.615
oldmarmot=rnorm(10^5)*sqrt(sigmas)+mus
hist(oldmarmot,nclass=200,col="gold",add=T)

# Monte Carlo CI

quantile(mus,c(.025,0.975))

# alternative method to get the CI from the exact formulas:
mean(maradult)+sd(maradult)*qt(0.975,31)/sqrt(32)
mean(maradult)-sd(maradult)*qt(0.975,31)/sqrt(32)

# getting a joint confidence region from simulations
toppoint=rep(0,10^5)
for (t in 1:10^5) toppoint[t]=topdensity(mus[t],sigmas[t])
boundpost=quantile(toppoint,.95)
plot(sigmas,mus,col="white")
points(sigmas[toppoint<boundpost],mus[toppoint<boundpost],pch=19,col="tomato")
points(sigmas[toppoint>boundpost],mus[toppoint>boundpost],pch=19,col="wheat")

#testing for mean equality on the weights of males and females
boxplot(marmotadult2$Peso~marmotadult2$Sex)
marmotadult2=marmot[marmot$Age=="ADULTO",]
marmotadult2m=marmotadult2[marmotadult2$Sex=="M",]
marmotadult2f=marmotadult2[marmotadult2$Sex=="F",]
wm=marmotadult2m$Peso
wf=marmotadult2f$Peso
zbar=sum(wm+wf)/32
zbar=(sum(wm)+sum(wf))/32
sz=sum((wm-zbar)^2)+sum((wf-zbar)^2)
sx=var(wm)*17
sy=var(wf)*13
(-sz+sx+sy-(zbar-3500)^2/(1/32+1/20)+(mean(wm)-3500)^2/(1/18+1/20)+(mean(wf)-3500)^2/(1/14+1/20))/(2*426^2)
-0.2807718+(log((1/32+1/20))-log((1/18+1/20))-log((1/14+1/20)))*.5-log(2*pi)*.5