######
###### hull
######
#!cp $SHOME/library/venables/morley.data .
system("cat morley.data")
#x11()
postscript(horiz=F)
x<-rnorm(50)
y<-rnorm(x)
hull<-chull(x,y)
plot(x,y)
polygon(x[hull],y[hull],den=15)
#hull.list<-list(x=x,y=y,hull=hull)
objects()
rm(x,y,hull)
######
###### regression
######
x<-1:20
w<-1+sqrt(x)/2
dummy<-data.frame(x=x,y=x+rnorm(x)*w)
dummy
fm<-lm(y~x, data=dummy)
summary(fm)
fm1<-lm(y~x,data=dummy,weight=1/w*2)
summary(fm1)
lrf<-loess(y~x,dummy)
attach(dummy)
plot(x,y)
lines(x,fitted(lrf))
abline(0,1,lty=3)
abline(coef(fm))
abline(coef(fm1),lty=4)
detach()
plot(fitted(fm),resid(fm),xlab="Fitted values", ylab="Residuals", main="Residuals vs Fitted")
qqnorm(resid(fm),main="Residuals Rankit Plot")
#regression.list<-list(fm=fm, fm1=fm1,lrf=lrf,dummy=dummy,w=w)
ls()
rm(x,w,fm,fm1,lrf,dummy)
######
###### morley
######
system("cat morley.data")
mm<-read.table("morley.data",header=T)
mm
attach(mm,1)
objects()
Expt<-factor(Expt)
Run<-factor(Run)
mm
detach()
attach(mm)
plot(Expt,Speed,main="Michaelson Morley Data",xlab="Experiment No.")
fm<-aov(Speed~Run+Expt,data=mm)
summary(fm)
fm0<-update(fm,.~.-Run)
summary(fm0)
anova(fm0,fm)
#morley.list<-list(mm=mm,fm=fm,fm0=fm0)
rm(mm,fm,fm0)
######
###### contour
######
x<-seq(-pi,pi,len=50)
y<-x
f<-outer(x,y,function(x,y)cos(y)/(1+x*2))
oldpar<-par()
#oldpar<-par(gr.state)
par(pty="s")
contour(x,y,f)
contour(x,y,f,nlevels=15,add=T)
fa<-(f-t(f))/2
contour(x,y,fa,nlevels=15)
par(oldpar)
#par(pty="m")
persp(x,y,f)
persp(x,y,fa)
#contour.list<-list(x=x,y=y,f=f,fa=fa,oldpar=oldpar)
rm(f,x,y,fa)
######
###### complex
######
th<-seq(-pi,pi,len=100)
z<-exp(1i*th)
par(pty="s")
plot(z,type="l")
w<-rnorm(100)+rnorm(100)*1i
w.wgt1<-w
w<-ifelse(Mod(w)>1,1/w,w)
plot(w,xlim=c(-1,1),ylim=c(-1,1),pch="+",xlab="x",ylab="y")
lines(z)
w.nonunif<-w
w<-sqrt(runif(100))*exp(2*pi*runif(100)*1i)
plot(w,xlim=c(-1,1),ylim=c(-1,1),pch="+",xlab="x",ylab="y")
lines(z)
w.unif<-w
plot(w.wgt1,xlim=c(-4,4),ylim=c(-4,4),pch="+",xlab="x",ylab="y")
points(w.unif,pch="!")
plot(w.nonunif,xlim=c(-1,1),ylim=c(-1,1),pch="+",xlab="x",ylab="y");lines(z)
plot(w.unif,xlim=c(-1,1),ylim=c(-1,1),pch="+",xlab="x",ylab="y");lines(z)
#complex.list<-list(w.wgt1=w.wgt1,w.unif=w.unif,w.nonunif=w.nonunif,th=th,z=z)
######
###### butterfly and exit
######
butterfly<-
function()
{
        theta <- seq(0, 24 * pi, len = 2000)
        radius <- exp(cos(theta)) - 2 * cos(4 * theta) + sin(theta/12)^5
        plot(radius * sin(theta), (-radius) * cos(theta), type = "l", 
                axes = F, xlab = "", ylab = "")
}
par(oldpar)
butterfly()
rm(oldpar)
#rm(z,th,w,w.wgt1,w.nonunif,w.unif)
graphics.off()
#q()