R Code 4: Integration + My first chain


1. Simulation is great.
1a. What is the variance of a uniform distribution?
Why not draw a whole bunch of samples and compute sample variance?
var(runif(10000))

1b. Let X1 and X2 be distributed as uniform(0,1). What does X1-X2 look like? What is the 25th percentile?
X1=runif(10000)
X2=runif(10000)
Y=X1-X2
hist(Y)
quantile(Y,p=.25)

#########################################
#Integration by Sampling

What is the area of a circle with radius 1?
my.col=c("red","blue")
M=100
x=runif(M,-1,1)
y=runif(M,-1,1)
#box is area 4.
distance.from.zero=sqrt(x*x+y*y)
in.circle=(distance.from.zero<1)
par(pty='s')
plot(x,y,pch=19,col=my.col[as.integer(in.circle)+1])
area=4*sum(in.circle)/M
#compare to pi*r^2

###################################################
dexg=function(t,mu,sigma,tau)
{
temp1=(mu/tau)+((sigma*sigma)/(2*tau*tau))-(t/tau)
temp2=((t-mu)-(sigma*sigma/tau))/sigma
(exp(temp1)*pnorm(temp2))/tau
}

mu=500
sigma=50
tau=200

rt=seq(0,1700,10)
plot(rt,dexg(rt,mu,sigma,tau),typ='l',lwd=2,col='darkblue')

#additive definition of ex-Gaussian
x ~ Exp(rate=1/tau)
y ~ Normal(mu,sigma^2)
z = x + z

#simuation method
M=20000
z=rexp(M,rate=1/tau)+rnorm(M,mu,sigma)
hist(z,prob=T,breaks=50)
lines(rt,dexg(rt,mu,sigma,tau))
####################################################

#conditional definition of ex-Gaussian
z|x ~ Normal(mu+x,sigma^2)
x ~ Exp(rate=1/tau)

#simulation from conditional form
x=rexp(M,rate=1/tau)
z=rnorm(M,x+mu,sigma)
hist(z,prob=T,breaks=50)
lines(rt,dexg(rt,mu,sigma,tau))

###################################################

#My first MCMC chain
library(MCMCpack)
#make data
dat=rnorm(4,110,10)
N=length(dat)
mean.dat=mean(dat)

#prior settings
mu.0=0
sig2.0=1000
a=1
b=1

M=1000
mu=1:M
sig2=1:M

#really bad starting values
mu[1]=0
sig2[1]=1

for (m in 2:M)
{
var.prime=1/((N/sig2[m-1])+(1/sig2.0))
mean.prime=((N*mean.dat/sig2[m-1])+(mu.0/sig2.0))*var.prime
mu[m]=rnorm(1,mean.prime,sqrt(var.prime))
s2.scale=sum((dat-mu[m])^2)/2+b
sig2[m]=rinvgamma(1,shape=a+N/2,scale=s2.scale)
}

#######################################################

plot(m
keep=10:M

hist(mu[keep])
hist(sig2[keep])

Research Methods 3020