postscript("hursh.ps")
par(cex.lab=2)
par(cex.axis=1.75)
par(mar=(c(5,4.7,4,2)+.1))
##dhursh<-read.csv("hursh.csv")
lmhursh<-lm(dhursh[,2]~dhursh[,1])
plot(dhursh[,1],dhursh[,2],main="",axes=FALSE, xlab="Diameter (microns)", ylab="Velocity (m/s)",xlim=c(0,20),ylim=c(0,120))
axis(1,pos=0)
axis(2,pos=0)
abline(as.numeric(lmhursh$coef[1]),as.numeric(lmhursh$coef[2]))
dev.off()


postscript("m1counts.ps")
par(cex.lab=2)
par(cex.axis=1.75)
par(mar=(c(5,4.7,4,2)+.1))
dm1counts<-read.csv("m1counts.csv",header=FALSE)
hist(dm1counts[,1],xlab="Spikes per Second",ylab="",main="",breaks=10,axes=FALSE,xlim=c(0,35),ylim=c(0,20))
axis(1,pos=0)
axis(2,pos=0)
dev.off()


postscript("rk2-ttest.ps")

par(mfrow=c(2,1))
drk=read.csv("rk.csv")
par(cex.lab=1.8)
par(cex.axis=1.6)

hist(drk[,1],breaks=10,xlim=c(0,30),ylim=c(0,6),xlab="SSSS",ylab="",axes=FALSE,main="")
axis(1,pos=0)
axis(2,pos=0)

hist(drk[,2],breaks=8,xlim=c(0,30),ylim=c(0,6),xlab="SSST",ylab="",axes=FALSE,main="")
axis(1,pos=0)
axis(2,pos=0)

dev.off()


postscript("hemi2.ps")
par(mfrow=c(2,1))
par(cex.lab=1.8)
par(cex.axis=1.6)

hist(hemi, breaks=40, xlim=c(0,1), xlab="Seconds", ylab = "", main=NULL, axes=FALSE)
axis(1, pos=0)
axis(2, pos=0)

hist(loghemi, breaks=40, xlim=c(-1.2,0), xlab="log(Seconds)", ylab = "", main=NULL, axes=FALSE)
axis(1, pos=0)
axis(2, pos=-1.2)


dev.off()



##exp-fill
lambda = 1
n = 1000

x.max = ceiling(qexp(.99, lambda))  # x grid displays at least 99% of area
x = seq(0, x.max, le=n)

fx = dexp(x, lambda)

## Shading options
pt = 2     # determine area to be shaded
np = 1000  # represents density of shaded area
shade.col = "lightgrey"

xp = seq(pt, max(x), le=np)
fxp = dexp(xp, lambda)

## Plotting starts here
postscript("exp-fill.ps")

par(cex.axis=1.75)
plot(x, fx, type="l", ylab="", xlab="", xlim=c(0,5), axes=F, lwd=2)
segments(xp,0,xp,fxp, col=shade.col)
axis(1, pos=0)
lines(x, fx)
segments(pt, 0, pt, dexp(pt, lambda))


## Plot labels
##title(expression(bold(P(X>2~"|"~lambda==1))))


## Saving Plot
##dev.copy2eps(file="exp-fill.ps")

dev.off()



postscript("qq.ps")
t3<-rt(200,df=3)
rn<-rnorm(200)
chi<-rchisq(200, 4)

par(mfrow=c(2,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
par(mar=(c(5,4.7,4,2)+.1))


qqnorm(rn, xlab="Theoretical Quantiles", ylab="", axes=F, main="", xlim=c(-3,3), ylim=c(-3,3))
axis(1,pos=-3)
axis(2,pos=-3)


qqnorm(chi, xlab="Theoretical Quantiles", ylab="", axes=F, main="")
axis(1,pos=0)
axis(2,pos=-3)

qqnorm(t3, xlab="Theoretical Quantiles", axes=F, ylab="", main="", ylim=c(-6,6))
axis(1,pos=-6)
axis(2,pos=-3)

dev.off()


postscript("qqnorm.ps")

par(mfrow=c(4,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
par(mar=(c(5,4.7,4,2)+.1))

for(i in 1:8)
{
norms<-rnorm(30)
qqnorm(norms, xlab="", ylab="", main="", axes = F, xlim=c(-3,3), ylim=c(-3,3)) 
axis(1,pos=-3)
axis(2,pos=-3)
}

dev.off()


postscript("for10.ps")

par(mfrow=c(1,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
par(mar=(c(5,4.7,4,2)+.1))

x5=c(0:4)
theta99=seq(.01,1,by=.01)

p5=x5
p99=theta99

for(i in 1:5)
{
p5[i]=dbinom(x5[i], 4, .5)
}

for(i in 1:99)
{
p99[i]=dbinom(1, 4, theta99[i])
}

plot(x5, p5, xlab="X", ylab="PDF", main="", axes=F, xlim=c(0,4), ylim=c(0,.4))
axis(1,pos=0)
axis(2,pos=0)
plot(theta99, p99, xlab=expression(theta), axes=F, ylab="PDF", main="", xlim=c(0,1), ylim=c(0,.5))
axis(1,pos=0)
axis(2,pos=0)
dev.off()


postscript("xbar-s.ps")

par(mfrow=c(1,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
par(mar=(c(5,4.7,4,2)+.1))

xbarData=matrix(nrow=1000, ncol=2)
for(i in 1:1000)
{
temp=rpois(100,10)
xbarData[i,1]=mean(temp)
xbarData[i,2]=var(temp)
}

hist(xbarData[,1], axes=F, main="", breaks=10, ylab="", xlab=expression(bar(Y)),xlim=c(4,16), ylim=c(0,300))
axis(1,pos=0)

hist(xbarData[,2], axes=F, main="", breaks=10, ylab="", xlab=expression(S^(2)),xlim=c(4,16), ylim=c(0,300))
axis(1,pos=0)

dev.off()

postscript("7resplot.ps")

d1=matrix(nrow=50, ncol=2)
for(i in 1:50)
{
d1[i, 2]=rnorm(1)
d1[i, 1]=runif(1, min=100, max=150)
}

d2=matrix(nrow=50, ncol=2)
for(i in 1:50)
{

d2[i, 1]=runif(1, min=100, max=150)
d2[i, 2]=.01*(d2[i,1]-125)^2-3+rnorm(1)
}



d3=matrix(nrow=50, ncol=2)
for(i in 1:50)
{
d3[i, 1]=runif(1, min=100, max=150)
d3[i, 2]=rnorm(1,mean=0, sd=(d3[i,1]*.008))

}

d4=matrix(nrow=50, ncol=2)
for(i in 1:50)
{
d4[i, 2]=rnorm(1)
d4[i, 1]=runif(1, min=100, max=150)
}
d4[50,2]=20



par(mfrow=c(2,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
par(mar=(c(5,4.7,4,2)+.1))

plot(d1[,1],d1[,2], axes=F, xlim=c(100,150), ylim=c(-3,3), main=NULL,xlab="", ylab="")
axis(1,pos=-3)
axis(2,pos=100)

plot(d2[,1],d2[,2], axes=F, xlim=c(100,150), ylim=c(-6,6),main=NULL,xlab="", ylab="")
axis(1,pos=-6)
axis(2,pos=100)

plot(d3[,1],d3[,2], axes=F, xlim=c(100,150), ylim=c(-4,4), main=NULL,xlab="", ylab="")
axis(1,pos=-4)
axis(2,pos=100)

plot(d4[,1],d4[,2], axes=F, xlim=c(100,150), ylim=c(-5,20), main=NULL,xlab="", ylab="")
axis(1,pos=-5)
axis(2,pos=100)

dev.off()



postscript("normbeta.ps")
par(cex.lab=2)
par(cex.axis=1.75)
par(mar=(c(5,4.7,4,2)+.1))

ppp=seq(.0005,.9995, by=.001)

##plot(cbind(ppp,pp), axes=F, xlab="",ylab="",cbind(dnorm(ppp,mean=.6,sd=.049),dbeta(pp,61,41)))


curve(dbeta(x,61,41),from=0,to=1,axes=F,xlab="",ylab="")
points(ppp,dnorm(ppp,mean=.6,sd=.049),type="p", cex=.1)
axis(1,pos=0)
dev.off()



postscript("QQp3359b.ps")

rn<-read.csv("p3359b.csv", header=F)


chi<-log10(rn)
t3<-(1/rn)
t4<-read.csv("p3437b.csv", header=F)

par(mfrow=c(2,2))
par(cex.lab=1.6)
par(cex.axis=1.4)
par(mar=(c(5,4.7,4,2)+.1))


qqnorm(rn[,1], xlab="Theoretical Quantiles", ylab="", axes=F, main="", xlim=c(-3,3), ylim=c(0,1))
axis(1,pos=0)
axis(2,pos=-3)


qqnorm(chi[,1], xlab="Theoretical Quantiles", ylab="", axes=F, main="", xlim=c(-3,3), ylim=c(-1.5,0))
axis(1,pos=-1.5)
axis(2,pos=-3)

qqnorm(t3[,1], xlab="Theoretical Quantiles", axes=F, ylab="", main="", ylim=c(0,15), xlim=c(-3,3))
axis(1,pos=-0)
axis(2,pos=-3)

qqnorm(t4[,1], xlab="Theoretical Quantiles", axes=F, ylab="", main="", ylim=c(0,1), xlim=c(-3,3))
axis(1,pos=-0)
axis(2,pos=-3)

dev.off()


source("flushplotold.r")
postscript("roc.ps")
par(lwd=2)
x<-seq(-4,6,length=200)
par(mfrow=c(2,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
matplot(x,cbind(dnorm(x,0,1),dnorm(x,2,1)),lwd=2,type="l",col="black",xlab="", ylab="",axes=FALSE)
axis(1,pos=0)


plot(1-pnorm(x,0,1),1-pnorm(x,2,1),type="l",lwd=2, xlab="", ylab="", axes=FALSE)
axis(1,pos=0)


matplot(x,cbind(dnorm(x,0,1),dnorm(x,1,1)),lwd=2,type="l",col="black", xlab="", ylab="", axes=FALSE)
axis(1,pos=0)



plot(1-pnorm(x,0,1),1-pnorm(x,1,1),type="l", lwd=2, xlab="", ylab="",axes=FALSE)
axis(1,pos=0)



dev.off()




## Author: Elan D. Cohen (edc@stat.cmu.edu)
## Date:   1/22/08
## Goal:   Simulate Galton father/son data (see Freedman, 78)


library(MASS)


## set bivariate normal parameters (from Freedman)
mu.x = 68
mu.y = 69
mu = c(mu.x, mu.y)
sd.x = 2.7
sd.y = 2.7
rho = 0.5
sig.xy = rho*sd.x*sd.y
Sig = matrix(c(sd.x^2, sig.xy, sig.xy, sd.y^2), nrow=2)

n = 1078
X = mvrnorm(n, mu, Sig)
X.lm = lm(X[,2]~X[,1])


## epsilon averaging
eps = 0.5
pt.1 = 64
mean.1 = mean(X[X[,1]>pt.1-eps & X[,1]<pt.1+eps,2])
pt.2 = 72
mean.2 = mean(X[X[,1]>pt.2-eps & X[,1]<pt.2+eps,2])


## Plotting starts here
postscript(file="galton.ps", width=9)
par(mar=(c(5,4.7,4,2)+.1))

plot(X[,2]~X[,1], xlab="Father's Height (inches)", ylab="Son's Height (inches)",
     pch=20, col="darkgray", cex=.5, cex.axis=1.75, cex.lab=2, font.lab=2)
abline(1,1,lty=2, lwd=2)  # draw y=x+1
abline(X.lm, lwd=2)  # draw regression line
abline(v=c(pt.1-eps, pt.1+eps, pt.2-eps, pt.2+eps), lty=2)  # draw vertical

points(pt.1, mean.1, pch="x", font=2, cex=1.5, col="black")
points(pt.2, mean.2, pch="x", font=2, cex=1.5, col="black")

dev.off()



### Galton-Contour Plot
f = function(x, y, mu=c(0,0), sd.x=1, sd.y=1, rho=0) {

  mu.x = mu[1]
  mu.y = mu[2]
  A = (x-mu.x)^2/sd.x^2 + (y-mu.y)^2/sd.y^2
  B = 2*rho/(sd.x*sd.y)*(x-mu.x)*(y-mu.y)
  return((A-B)/(1-rho^2))
}


N = 100
x = seq(min(X[,1]), 1.005*max(X[,1]), le=N)
y = seq(min(X[,2]), max(X[,2]), le=N)
z = outer(x, y, f, mu=mu, sd.x=sd.x, sd.y=sd.y, rho=rho)

p.seq = c(.8, .9, .95, .99)
cont.lev = qchisq(p.seq, 2)
cont.lab = c("80%", "90%", "95%", "99%")


## Plotting starts here
postscript(file="galton-contour.ps", width=9)

par(mar=(c(5,4.7,4,2)+.1))
contour(x, y, z,levels=cont.lev, labcex=1, labels=cont.lab,
        xlab="Father's Height (inches)", ylab="Son's Height (inches)",
        cex.axis=1.75, font.lab=2, cex.lab=2, ylim=c(56,81), xlim=c(58,77))

abline(1,1,lty=2, lwd=2)  # draw y=x+1
abline(v=c(pt.1-eps, pt.1+eps, pt.2-eps, pt.2+eps), lty=2)  # draw vertical
## True Regresion Line: E(Y|X)
b = rho*sd.y/sd.x  # slope
a = mu.y - mu.x*b  # intercept
abline(a, b, lwd=2)

dev.off()




### EOF ###





## Author:  Elan Cohen (edc@stat.cmu.edu)
## Date:    1/22/08
## Goal:    Produce example perspective and contour plots of the bivariate
##          Normal distribution.
## Run:     > source("bivariate-normal.r")


### Start: Function definitions ###

## Calculated bivariate normal density at /one/ point (x,y).
## Must pass full covariance matrix (Sig), or sd1, sd2 and rho.
## See <dmvnorm> below.
fxy = function(x, y, mu, Sig, sd1, sd2, rho) {
  
  if(missing(mu)) mu=c(0,0)
  
  if(!missing(Sig)) {
    sd1 = sqrt(Sig[1,1])
    sd2 = sqrt(Sig[2,2])
    if(Sig[1,2] != Sig[2,1]) {
      print("Covariance matrix is not symmetric... Returning <NULL>.")
      return(NULL)
    }
    rho = Sig[1,2]/(sd1*sd2)
  }
  else if(missing(rho) || missing(sd1) || missing(sd2)) {
    sd1 = sd2 = 1
    rho = 0
  }

  Q = (x-mu[1])^2/sd1^2 + (y-mu[2])^2/sd2^2 -
    2*rho*(x-mu[1])*(y-mu[2])/(sd1*sd2)
  
  1/(2*pi*sd1*sd2*sqrt(1-rho^2))*exp(-Q/(2*(1-rho^2)))
}


## Calls persp() with preferred arguments
persp.plot = function(x, y, z, main="Bivariate Normal Density",
  theta=30, phi=25, r=50, d=.1, expand=0.5, ltheta=90, lphi=180,
  shade=0.5, ticktype="simple", nticks=5, col="black", zlab="", ...) {
  
  persp(x, y, z, main=main,
        theta=theta, phi=phi, r=r, d=d, expand=expand, ltheta=ltheta,
        lphi=lphi, shade=shade, ticktype=ticktype, nticks=nticks,
        col=col, zlab=zlab, ...)
}


## Creates covariance matrix from sd.x, sd.y, and rho
calc.Sig = function(sd.x, sd.y, rho) {

  sig.xy = rho*sd.x*sd.y
  matrix(c(sd.x^2, sig.xy, sig.xy, sd.y^2), nrow=2)
}


## Returns bivariate normal density for specified x-y grid
dmvnorm = function(x, y, mu, Sig) {

  if(missing(mu)) mu = c(0,0)
  if(missing(Sig)) Sig = diag(2)

  outer(x, y, fxy, mu, Sig)
}


## This is only the kernel of the bivariate Normal density
## x is a 2x1 vector
f = function(x, y, mu=c(0,0), sd.x=1, sd.y=1, rho=0) {

  #t(X-mu)%*%solve(Sig)%*%(X-mu)
  mu.x = mu[1]
  mu.y = mu[2]
  A = (x-mu.x)^2/sd.x^2 + (y-mu.y)^2/sd.y^2
  B = 2*rho/(sd.x*sd.y)*(x-mu.x)*(y-mu.y)
  return((A-B)/(1-rho^2))
}


### End: Function definitions ###



### Create perspective and contour plots

## The value of N affects the density of lines and hence the darkness.
## Unfortunately, small values of N result in a rough plot, while large values
## result in a dark plot.  This is also dependent on the size of the X window
## (a large window size will appear lighter than a smaller window size for a
## fixed N).  If 'border=NA' is set, these lines don't appear.
N = 100
x = y = seq(-3.2,3.2,le=N)  # create x-y grid of size NxN
mu = c(0,0)


## Define sequence of parameters
n.plot = 6
rho.seq = rep(c(0, 0.75, -0.75), 2)
sd.x.seq = rep(c(1, .5), each=n.plot/2)
sd.y.seq = rep(1, n.plot)

expr.seq = c(expression(list(sigma[x]==sigma[y], ~rho==0)),
  expression(list(sigma[x]==sigma[y], ~rho==0.75)),
  expression(list(sigma[x]==sigma[y], ~rho==-0.75)),
  expression(list(2*sigma[x]==sigma[y], ~rho==0)),
  expression(list(2*sigma[x]==sigma[y], ~rho==0.75)),
  expression(list(2*sigma[x]==sigma[y], ~rho==-0.75)))

p.seq = c(.8, .9, .95, .99)
cont.lev = qchisq(p.seq, 2)
cont.lab = c("80%", "90%", "95%", "99%")


## Plotting starts here
postscript(file="bivariate-normal-%01d.ps", horizontal=F,
           width=7, height=10.5, onefile=F)  # 1.5:1 aspect ratio


par(mfrow=c(3,2))  # 1.5:1 aspect ratio
par(cex.axis=1.55)
par(cex.lab=1.75)

for(j in 1:n.plot) {

  ## Perspective Plot
  z = dmvnorm(x, y, mu, calc.Sig(sd.x.seq[j], sd.y.seq[j], rho.seq[j]))
  persp.plot(x, y, z, main="", col="lightgrey", border=NA,
             axes=F)
  mtext(expr.seq[j],cex=1.4)

  ## Contour Plot
  z = outer(x, y, f, mu, sd.x.seq[j], sd.y.seq[j], rho.seq[j])
  contour(x, y, z, levels=cont.lev, labels=cont.lab,
          xlab="x", ylab="y")
  abline(v=0, h=0, lty=3, col="darkgrey")

}

dev.off()




### EOF ###





source("flushplot.r") 

postscript("4dens.ps")
par(mfrow=c(2,2))
par(cex.lab=1.8)
par(cex.axis=1.6)
par(lwd=2)

xdens<- c(seq(.00001,8,length.out=500))
ydens<- dexp(xdens)
flushplotold.fun(xdens,ydens,xlab="")

xdens<- c(seq(.00001,10,length.out=500))
ydens<- dgamma(xdens,2)
flushplotold.fun(xdens,ydens,xlab="")

xdens<- c(seq(.00001,1,length.out=500))
ydens<- dbeta(xdens,6,6)
flushplotold.fun(xdens,ydens,xlab="")

xdens<- c(seq(-5,5,length.out=500))
ydens<- dnorm(xdens)
flushplotold.fun(xdens,ydens,xlab="")


dev.off()



source("flushplot.r") 

postscript("invgauss.ps")
%par(mfrow=c(2,2))
par(lwd=2)
par(cex.lab=1.8)
par(cex.axis=1.6)
dinvgauss<-function(x,mu,lambda){
  sqrt(lambda/(2*pi*x^3))*exp(-lambda*(x-mu)^2/(2*x*mu^2))
}
xdens<- c(seq(.00001,10,length.out=500))
ydens<- dgamma(xdens,2)
y2dens<-dinvgauss(xdens,2,8)
flushplotold.fun(xdens,y2dens,xlab="",lwd=2,lty="dashed")
lines(xdens,ydens,lwd=2)

dev.off()



postscript("hemi.ps")
par(mfrow=c(2,2))

par(font=1)
par(cex.axis=1.6)
par(cex.lab=1.8)
hist(hemi,breaks=c(0,seq(0:10)/10),main=NULL,axes=FALSE,xlim=c(0,1),xlab="Time (s)",ylab=NULL, lwd=2)
axis(1,pos=0)
axis(2,pos=0)
hist(hemi,breaks=c(0,seq(0:50)/50),main=NULL,axes=FALSE, xlim=c(0,1), xlab="Time (s)",ylab=NULL, lwd=2)
axis(1,pos=0)
axis(2,pos=0)
hist(hemi,breaks=c(0,seq(0:20)/20),main=NULL,axes=FALSE, xlim=c(0,1), xlab="Time (s)",ylab=NULL, lwd=2)
axis(1,pos=0)
axis(2,pos=0)
hist(hemi,breaks=c(0,seq(0:20)/20+.025),main=NULL,axes=FALSE, xlim=c(0,1), cex.lab=1.5, xlab="Time (s)",ylab=NULL, lwd=2)
axis(1,pos=0)
axis(2,pos=0)
dev.off()



postscript("normbinom.ps")
par(cex.lab=2)
par(cex.axis=1.75)
xvals=c(10:70)
plot(xvals,dbinom(xvals,p=.4,size=100),xlab="",ylab="",lwd=2,axes=FALSE)
lines(xvals,dnorm(xvals,mean=40,sd=sqrt(24)),col="black",lwd=2)
axis(1,pos=0)
dev.off()




      






source("flushplotOLD.r") 


postscript("binom.ps")
par(mfrow=c(2,2))
par(cex.lab=2)
par(cex.axis=1.75)
par(cex.main=2)
par(cex.lab=2)

xdens<- c(seq(0,4,length.out=5))
ydens<- dbinom(xdens,4,.4)
plot(xdens/4,ydens,axes=FALSE,xlab="",ylab="",main="")
axis(1,pos=0)

xdens<- c(seq(0,10,length.out=11))
ydens<- dbinom(xdens,10,.4)
plot(xdens/10,ydens,axes=FALSE,xlab="",ylab="",main="")
axis(1,pos=0)

xdens<- c(seq(0,25,length.out=26))
ydens<- dbinom(xdens,25,.4)
plot(xdens/25,ydens,axes=FALSE,xlab="",ylab="",main="")
axis(1,pos=0)

xdens<- c(seq(0,100,length.out=101))
ydens<- dbinom(xdens,100,.4)
plot(xdens/100,ydens,axes=FALSE,xlab="",ylab="",main="")
axis(1,pos=0)

dev.off()




source("flushplotOLD.r") 


postscript("chisq.ps")
par(mfrow=c(2,2))
par(cex.lab=2)
par(cex.axis=1.75)
xdens<- c(seq(.00001,.5,length.out=500))
ydens<- dchisq(xdens,1)
flushplot.fun(xdens,ydens,xlab="")

xdens<- c(seq(.00001,12,length.out=500))
ydens<- dchisq(xdens,4)
flushplot.fun(xdens,ydens,xlab="")

xdens<- c(seq(.00001,25,length.out=500))
ydens<- dchisq(xdens,9)
flushplot.fun(xdens,ydens,xlab="")

xdens<- c(seq(.00001,35,length.out=500))
ydens<- dchisq(xdens,16)
flushplot.fun(xdens,ydens,xlab="")

dev.off()