# plotting a bunch of gamma densities just to see what they look like
# source one line at the time

xx _ seq(.01,9.99,length=100)
plot(xx, dgamma(xx,1), lwd=3, type="l")
lines(xx, dgamma(xx,2), lwd=3, type="l", col=2)
lines(xx, dgamma(xx,4), lwd=3, type="l", col=3)
lines(xx, .5*dgamma(xx*.5,2), lwd=3, type="l", col=4)
lines(xx, 8*dgamma(xx*8,16), lwd=3, type="l", col=5)

# gamma updating
# assumes nu=1, a0=16, a=6, 

plot(c(0,10),c(0,.5),type="n",xlab="THETA",ylab="DENSITY")
lines(xx, dgamma(xx,6,1), lwd=3)
points(5,0,cex=2); points(4.2,0,cex=2,col=2)
lines(xx, dgamma(xx,6+16,1/(4.2+1)), lwd=3, col=2)

# simulate points from two independent gamma and add some normal
# measurement error  

plot(c(0,4),c(0,4),type="n",xlab="GREEN",ylab="RED") # sets stage
abline(0,1)                                          # x==y line
x_rgamma(1,2); y_rgamma(1,2); points(x,y, pch=1, cex=2); # one black point
xx _ rnorm(10,x,.1); yy _ rnorm(10,y,.1); points(xx,yy,col=4) # blue pts aroun it
for (i in 1:10) abline(0,yy[i]/xx[i], col=4)                  # ratios
xx _ rnorm(10,x,.1); yy _ rnorm(10,y,.1); points(xx,yy,col=6)
for (i in 1:10) abline(0,yy[i]/xx[i], col=6)