R Stuff, Chapter 2


1. Functions

#order vs labels
y=rnorm(10)
y=rnorm(10,10)
y=rnorm(10,0,3)
y=rnorm(sd=3,mean=10,n=40)
#often, you will need labels

#local vs. global
x=1:3
y=4:6

widget=function(w,z) sum(w^2+z)

widget(w=x,z=y)
widget(z=y,w=x)
widget(x,y)

widget=function(w,z)
{
x=sum(w^2+z) #this x is different than 1:3
return(x-2)
}

w=x
z=y

widget(w=w,z=z)
widget(z,w)
widget(w=z,z=w)

2. Logs in Plots

year=c(1820,1860,1880,1920,1940,1980,2000)
pirates=c(35000,45000,20000,15000,5000,400,17)
temp=c(14.2,14.4,14.5,14.7,15.0,15.4,15.8)

#without logs
plot(pirates,temp,xlim=c(0,50000),
xlab="Number of Pirates",ylab="Global Temperature (celsius)")
text(x=pirates,y=temp,labels=year,adj=c(-.2,.5))

#with logs
plot(log10(pirates),temp,xlim=c(1,5),
xlab="Log of Number of Pirates (base 10)",ylab="Global Temperature (celsius)")
text(x=log10(pirates),y=temp,labels=year,adj=c(-.2,.5))

#prettier
plot(log10(pirates),temp,axes=F,xlim=c(1,5),typ='n',
xlab="Number of Pirates",ylab="Global Temperature (celsius)")
points(log10(pirates),temp,cex=1.3,pch=21,bg='cyan')
axis(2)
axis(1,at=1:5,labels=c(10,100,1000,10000,100000))
box()
text(x=log10(pirates),y=temp,labels=year,adj=c(-.2,.5))

#My Preference
plot(log10(pirates),temp,axes=F,xlim=c(1,5),typ='n',
xlab="Number of Pirates",ylab="Global Temperature (celsius)")
axis(2)
axis(1,at=1:5,labels=c(10,100,1000,10000,100000))
box()
text(x=log10(pirates),y=temp,labels=year,adj=c(.5,.5))

3. Likelihood of Binomial
likelihood=function(p,y,N) dbinom(y,N,p)

#Look at likelihood
p=seq(0,1,.01)
like=likelihood(p,5,10)
plot(p,like,type='l') #why typ='l' and not typ='h'

#fancy code for Jeff
p=seq(0,1,.01)
y=0:10
f=function(y,p,N) dbinom(y,N,p)
a=outer(y,p,f,N=10)
image(y,p,a,col = gray((0:128)/128))
image(y,p,log(a),col = gray((0:128)/128))

#maximize likelihood
loglike=function(p,y,N)
return(dbinom(y,N,p,log=T))

optimize(loglike,interval=c(0,1),maximum=T,y=y,N=N)

#Case 2, is the coin fair?
N=20
y=13
mle.general=y/N
mle.restricted= .5 #by definition
log.like.general=dbinom(y,N,mle.general,log=T)
log.like.restricted=dbinom(y,N,mle.restricted,log=T)

4. Simulation Method For G^2

I have asserted that G^2 ~ Chi-Square(1) for the fair toast problem. I show it using the simulation method.

So, lets assume the null is true:

N=20
y=rbinom(100000,N,.5)
log.like.general=dbinom(y,N,y/N,log=T)
log.like.restricted=dbinom(y,N,.5,log=T)
g2=2*(log.like.general-log.like.restricted)
hist(g2)
quantile(g2,.95)

Why so bad; N is too small and everything is too discretized. Let's try N-200

N=200
y=rbinom(100000,N,.5)
log.like.general=dbinom(y,N,y/N,log=T)
log.like.restricted=dbinom(y,N,.5,log=T)
g2=2*(log.like.general-log.like.restricted)
hist(g2)
quantile(g2,.95)

AttachmentSize
PiratesVsTemp Original43.51 KB
DeathsInBible Original16.27 KB

Research Methods 3020