CTFSworkshop()
CTFSworkshop(data=TRUE)
% system("~/meetings_workshops/Rmodeling/CreateCourse.bat outline Tupper2012 Public/Workshops")
system("~/meetings_workshops/Rmodeling/CreateCourse.bat outline Tupper2012 Public/Workshops")
system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops")
system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops")
system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops")
search()
ls(2)
dim(2)
dim(agb)
head(agb)
plot(agb~ht,data=agb)
plot(I(log(agb))~I(log(ht)),data=agb)
lm(agb~ForestType+ht,data=agb)
ls(2)
ditch(2)
attach('~/Desktop/agb.rdata')
ls(2)
head(agb)
lm(agb~ForestType+ht,data=agb)
lm(agb~ForestType+ht,data=agb)
agb$logdbh=log(agb$dbh)
agb$logagb=log(agb$agb)
TreeMass=agb
save(TreeMass,file='~/Desktop/TreeMass.rdata')
head(TreeMass)
table(TreeMass$ForestType)
tapply(TreeMass$Rainfall,TreeMass$ForestType,range)
tapply(TreeMass$Alt,TreeMass$ForestType,range)
weird=TreeMass$ForestType=='Moist-'
table(weird)
TreeMass$ForestType[weird]
weird=TreeMass$ForestType=='Moist-' & !is.na(TreeMass$ForestType
)
weird=TreeMass$ForestType=='Moist-' & !is.na(TreeMass$ForestType)
TreeMass$ForestType[weird]
TreeMass$ForestType[weird]='Moist'
save(TreeMass,'~/meetings_workshops/Rmodeling/treemass.rdata')
save(TreeMass,'~/meetings_workshops/Rmodeling/treemass.rdata')
treemass=TreeMass
save(treemass,file='~/meetings_workshops/Rmodeling/treemass.rdata')
head(treemass)
lm(logagb~logdbh+ForestType,data=treemass)
fit=lm(logagb~logdbh+ForestType,data=treemass)
coef(fit)
abline(coef(fit)[1:2])
abline(coef(fit)[1],sum(coef(fit)[2:3]))
abline(sum(coef(fit)[c(1,3)]),coef(fit)[2])
ls(2)
plot(logagb~logdbh,data=treemass)
abline(sum(coef(fit)[c(1,3)]),coef(fit)[2])
abline(sum(coef(fit)[c(1,3)]),coef(fit)[2],col='red')
sum(coef(fit)[c(1,3)])
abline(coef(fit)[1:2],col='red')
abline(sum(coef(fit)[c(1,3)]),coef(fit)[2],col='green')
abline(sum(coef(fit)[c(1,4)]),coef(fit)[2],col='blue')
plot(logagb~logdbh,data=treemass,cex=.5,pch=16)
abline(sum(coef(fit)[c(1,4)]),coef(fit)[2],col='blue')
abline(sum(coef(fit)[c(1,3)]),coef(fit)[2],col='green')
abline(coef(fit)[1:2],col='red')
treemass$ctrx=treemass$logdbh-mean(treemass$logdbh)
plot(logagb~ctrx,data=treemass,cex=.5,pch=16)
fit=lm(logagb~ctrx+ForestType,data=treemass)
abline(coef(fit)[1:2],col='red')
coef(fit)
system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops")
system("~/meetings_workshops/Rmodeling/CreateCourse.bat outline Tupper2012 Public/Workshops")
treemass$ctrx=treemass$logdbh-mean(treemass$logdbh)
plot(logagb~logdbh,data=treemass,cex=.5,pch=16)
fit=lm(logagb~logdbh,data=treemass)
coef(fit)
abline(coef(fit))
abline(coef(fit),col='red')
fit=lm(logagb~ctrx+ForestType,data=treemass)
abline(fit)
fit=lm(logagb~logdbh,data=treemass)
plot(logagb~logdbh,data=treemass,cex=.5,pch=16)
fit=lm(logagb~logdbh,data=treemass)
search()
options(width=125)
search()
options(width=80)
search()
ls(6)
ls()
ls(6)
search()
ls(6)
head(treeht)
search()
ls(2)
head(treeheight)
head(treemass)
search()[1:4]
attach("~/meetings_workshops/Rmodeling/treemass.rdata")
rm(treemass)
ditch(2)
attach("~/meetings_workshops/Rmodeling/treemass.rdata")
ls(2)
head(treemass)
options(width=60)
head(treemass)
options(width=65)
head(treemass)
ls(2)
detach(2)
attach("~/meetings_workshops/Rmodeling/treemass.rdata")
ls(2)
head(treemass)
fit=lm(logagb~logdbh,data=treemass)
plot(logagb~logdbh,data=treemass,cex=.5,pch=16)
plot(treemass$logagb~treemass$logdbh,cex=.5,pch=16)
summary(fit)
coef(fit)
fit=lm(logagb~logdbh,data=treemass)
coef(fit)
abline(fit)
abline(fit,col='red')
plot(logagb~logdbh,data=treemass,cex=.5,pch=16)
fit=lm(logagb~logdbh,data=treemass)
abline(fit)
abline(fit,col='red')
coef(fit)
source('~/meetings_workshops/Rmodeling/teaching.functions.r')
args(sumsq.linearmodel)
sumsq.linearmodel(beta=coef(fit),x=treemass$logdbh,y=treemass$logagb)
sumsq.linearmodel(beta=coef(fit),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(5,0),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(5,2),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(5,2),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(-5,2),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(0,.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(0,.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
sumsq.linearmodel(beta=c(0,1.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
?optim
?optim
optim(par=c(0,0),fn=sumsq.linearmodel)
?optim
optim(par=c(0,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb))
optim(par=c(0,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb)
sumsq.linearmodel(beta=c(0,1.5),x=treemass$logdbh,y=treemass$logagb,graphit=TRUE)
?optim
optim(par=c(0,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb)
coef(fit)
optim(par=c(10,0),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb)
optim(par=c(10,-10),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb)
?optim
?optim
control=list(fnscale=(-1))q
?optim
llike.linearmodel(param=c(-2,2),x=treeheight$logdbh,y=treeheight$logdbh)
dnorm(1,mean=0,sd=1)
dnorm(0:5,mean=0,sd=1)
dnorm(-5:5,mean=0,sd=1)
plot(-5:5,dnorm(-5:5,mean=0,sd=2)
)
plot(-5:5,dnorm(-5:5,mean=0,sd=2,log=TRUE))
?dnorm
plot(-5:5,dnorm(-5:5,mean=1:11,sd=2,log=TRUE))
dnorm(-5:5,mean=1:11,sd=2,log=TRUE)
dnorm(-5:5,mean=1,sd=2,log=TRUE)
dnorm(-5:5,mean=1:11,sd=1:11,log=TRUE)
llike.linearmodel(param=c(-2,2),x=treeheight$logdbh,y=treeheight$logdbh)
llike.linearmodel(param=c(-2,2,1),x=treeheight$logdbh,y=treeheight$logdbh)
llike.linearmodel(param=c(-2,2,1),x=treemass$logdbh,y=treemass$logdbh)
debug(llike.linearmodel)
llike.linearmodel(param=c(-2,2,1),x=treemass$logdbh,y=treemass$logdbh)
param
summary(x)
SD
summary(yhat)
plot(x,y)
Q
llike.linearmodel(param=c(-2,2,1),x=treemass$logdbh,y=treemass$logagb)
plot(x,y)
points(x,yhat,col='red')
length(llike)
plot(x,llike)
x11()
points(x,yhat,col='red')
plot(x,y)
points(x,yhat,col='red')
plot(y,yhat,col='red')
plot(y,llike,col='red')
total
optim(par=c(10,-10),fn=sumsq.linearmodel,x=treemass$logdbh,y=treemass$logagb)
optim(par=c(10,-10),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb)
Q
undebug(llike.linearmodel)
optim(par=c(10,-10),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb)
optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb)
optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1)))
coef(fit)
llike.linearmodel(param=c(-2.14,2.51,.43),x=treemass$logdbh,y=treemass$logagb)
llike.linearmodel(param=c(-2.14,2.51,1),x=treemass$logdbh,y=treemass$logagb)
llike.linearmodel(param=c(-2.14,2.51,.1),x=treemass$logdbh,y=treemass$logagb)
llike.linearmodel(param=c(-2.14,2.51,.42),x=treemass$logdbh,y=treemass$logagb)
llike.linearmodel(param=c(-2.14,2.51,.41),x=treemass$logdbh,y=treemass$logagb)
llike.linearmodel(param=c(-2.14,2.51,.44),x=treemass$logdbh,y=treemass$logagb)
optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1)))
plot(logagb~logdbh,data=treemass,cex=.5,pch=16)
fit=lm(logagb~logdbh,data=treemass)
sd(resid(fit))
abline(fit,col='red')
summary(fit)
optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1)))
?optim
search()
ls(5)
head(cecrin)
search()
opt=optim(par=c(10,-10,1),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb,control=list(fnscale=(-1)),method='SANN')
opt$par
opt$par
optim(par=c(10,-10),fn=llike.linearmodel,x=treemass$logdbh,y=treemass$logagb)
fit1=optim(par=c(0,0,0),fn=llike.linearmodel,x=cerin$dbh1,y=cecrin$gr12)
fit1=optim(par=c(0,0,0),fn=llike.linearmodel,x=cerin$dbh1,y=cecrin$gr12)
fit1=optim(par=c(0,0,0),fn=llike.linearmodel.full,x=cerin$dbh1,y=cecrin$gr12)
fit1=optim(par=c(0,0,0.1),fn=llike.linearmodel.full,x=cerin$dbh1,y=cecrin$gr12)
fit1=optim(par=c(0,0,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12)
fit1$par
fit1=optim(par=c(0,0,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1)))
fit1$par
fit1$value
plot(gr12~dbh1,data=cecrin)
abline(fit1$par[1:2])
fit1=optim(par=c(0,-1,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1)))
fit1$par
abline(fit1$par[1:2])
fit1=optim(par=c(0,-1,0.1),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1)),method='SANN')
fit1$par
abline(fit1$par[1:2])
fit1$par
fit1$value
fit2=optim(par=fit1$par,fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1)))
fit2$par
fit2$value
coef(lm(gr12~dbh1,data=cecrin))
abline(fit2$par[1:2],'red')
abline(fit2$par[1:2],col='red')
fit3=optim(par=c(fit2$par,0),fn=llike.linearmodel.full,x=cecrin$dbh1,y=cecrin$gr12,control=list(fnscale=(-1)))
fit3$value
abline(fit3$par[1:2],col='blue')
fit3$par
abline(fit3$par[3:4],col='blue',lty='dashed')
abline(sum(fit3$par[1:2],fit3$par[3:4]),col='blue')
fit3$par[1:2]+fit3$par[3:4])
fit3$par[1:2]+fit3$par[3:4]
fit3$par
abline(fit3$par[1:2]+fit3$par[3:4],col='blue')
plot(gr12~dbh1,data=cecrin)
abline(fit3$par[1:2],col='blue')
abline(fit3$par[1:2]+fit3$par[3:4],col='blue')
abline(fit3$par[1:2]+2*fit3$par[3:4],col='blue')
abline(fit3$par[1:2]-2*fit3$par[3:4],col='blue')
fit2$par
plot(gr12~dbh1,data=cecrin)
abline(fit2$par[1:2],col='blue')
abline(c(sum(fit2$par[c(1,3)]),fit2$par[2]),col='blue')
abline(fit3$par[1:2]+2*fit3$par[3:4],col='blue')
abline(fit3$par[1:2]+fit3$par[3:4],col='blue')
system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops")
source('~/meetings_workshops/Rmodeling/RfuncBayesWorkshop.r')
source('~/meetings_workshops/Rmodeling/RcmdBayesWorkshop.r')
Q
graphics.off()
sample=10
weights=rnorm(n=sample,mean=mu,sd=sigma)
# If the mean were 125 and the SD 25, the likelihood is
dnorm(x=weights,mean=mu,sd=sigma)
# There are 10 points so 10 likelihoods, and it is easier to work with log-likelihood from the start.
# We always need the sum of the log-likelihoods, which is the (log) likelihood of the entire set of data. The absolute
# probability is not useful. 
loglike=dnorm(x=weights,mean=mu,sd=sigma,log=TRUE)
totalLogLike=sum(loglike)
cat(1)
cat('1','\n')
theta=0.76
N=10
dbinom(x=10,size=N,prob=theta)
.76^10
dbinom(x=9,size=N,prob=theta)
dbinom(x=0:10,size=N,prob=theta)
dbinom(x=11,size=N,prob=theta)
like=dbinom(x=11,size=N,prob=theta)
sum(like)
like=dbinom(x=0:10,size=N,prob=theta)
sum(like)
S=8
N
dbinom(x=S,size=N,prob=.2)
dbinom(x=S,size=N,prob=.3)
dbinom(x=S,size=N,prob=.4)
dbinom(x=S,size=N,prob=.9)
dbinom(x=S,size=N,prob=.8)
llike.binom=function(par,S,N) return( dbinom(x=S,size=N,prob=par,log=TRUE))
llike.binom
?optimize
optimize(f=llike.binom,N=N,S=S)
optimize(f=llike.binom,N=N,S=S,interval=c(0,1))
optimize(f=llike.binom,N=N,S=S,interval=c(0,1),maximum=TRUE)
optimize(f=llike.binom,N=10,S=8,interval=c(0,1),maximum=TRUE)
optimize(f=llike.binom,N=10,S=7,interval=c(0,1),maximum=TRUE)
7/10
theta=seq(0,1,by=0.01)
postlike=dbinom(x=S,size=N,prob=theta)
dev.set(2)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
theta
S
N
postlike=dbinom(x=80,size=100,prob=theta)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
x11()
postlike=dbinom(x=8,size=10,prob=theta)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
optimize(f=llike.binom,N=10,S=8,interval=c(0,1),maximum=TRUE)
optimize(f=llike.binom,N=100,S=80,interval=c(0,1),maximum=TRUE)
optimize(f=llike.binom,N=1000,S=800,interval=c(0,1),maximum=TRUE)
postlike=dbinom(x=8,size=10,prob=theta)
postlike=dbinom(x=800,size=1000,prob=theta)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
postlike=dbinom(x=8,size=10,prob=theta)
plot(theta,postlike,pch=16,type="l",main="binomial posterior",ylab="prob")
theta=0.76
dbinom(x=S,size=N,prob=theta)
theta=0.74
dbinom(x=S,size=N,prob=theta)
theta=0.72
dbinom(x=S,size=N,prob=theta)
theta=0.78
dbinom(x=S,size=N,prob=theta)
theta=0.82
dbinom(x=S,size=N,prob=theta)
S
N
dbinom(x=S,size=N,prob=.8)
postlike=dbinom(x=S,size=N,prob=theta)
sum(postlike)
k=sum(0.01*postlike)
post=postlike/k
sum(0.01*post)
lines(theta,post,type="l",col="red",main="adjusted binomial posterior",ylab="prob")
dbinom(8,10,prob=.19)
dbinom(8,10,prob=.20)
dbinom(8,10,prob=.01)
showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=100)
showBinomialMCMC(thetastart=.2,S=8,N=10)
Q
showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=100)
graphics.off()
showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=1000)
thetaMCMC=showBinomialMCMC(thetastart=.2,S=8,N=10,pause=FALSE,steps=1000,showstep=4000)
summary(thetaMCMC)
quantile(thetaMCMC,prob=.95)
quantile(thetaMCMC,prob=.05)
abline(v=.931,col='red')
abline(v=.494,col='red')
dev.set(3)
abline(v=.931,col='red')
abline(v=.494,col='red')
sd(thetaMCMC)
var(thetaMCMC)
median(thetaMCMC)
mean(thetaMCMC)
length(thetaMCMC)
length(whichMCMC<.6)
length(which(thetaMCMC<.6))
graphics.off()
fit3$value
system("~/meetings_workshops/Rmodeling/CreateCourse.bat assignments Tupper2012 Public/Workshops")
save.history(file='~/meetings_workshops/Rmodeling/Tupper2012/Tues7May.txt')
savehistory(file='~/meetings_workshops/Rmodeling/Tupper2012/Tues7May.txt')