library(MCMCpack)
#data(women)
#names(women) <- c("h","w")
#generate fake data: y = 1+2x-3x^2 + e
n <- 200
e <- rnorm(n)
x <- runif(n,0,2)
y <- 1+2*x-3*x^2+e

plot(x,y)
mydata <- data.frame(x=x,y=y)

library(MCMCpack)
# params for MCMCregress:
cc0=1; dd0=1000; BB0=1e-5
# straight line with intercept
m1 <- MCMCregress(y~x,data=mydata,c0=cc0,d0=dd0,B0=BB0,marginal.likelihood="Chib95")
sm1 <- summary(m1)
abline(sm1$stat[1:2,1],col="red")
# horiz. line
m0 <- MCMCregress(y~1,data=mydata,c0=cc0,d0=dd0,B0=BB0,marginal.likelihood="Chib95")
sm0 <- summary(m0)
abline(h=sm0$stat[1,1],col="blue")

# quadratic
m2 <- MCMCregress(y~x+I(x^2),data=mydata,c0=cc0,d0=dd0,B0=BB0,marginal.likelihood="Chib95")
sm2 <- summary(m2)
p2 <- function(x) {
    a <- sm2$stat[1:3,1]
    a[1]+a[2]*x+a[3]*x^2
}
curve(p2(x),add=TRUE)

# Bayes Factors for Model Selection
BF <- BayesFactor(m0,m1,m2)
PPM <- PostProbMod(BF)

PPM




Bregs <- function(cc0=1,dd0=1000,BB0=1e-5) {
	# no intercept
	m0 <- MCMCregress(y~0+x,data=mydata,c0=cc0,d0=dd0,B0=BB0,marginal.likelihood="Chib95")
	# with intercept
	m1 <- MCMCregress(y~x,data=mydata,c0=cc0,d0=dd0,B0=BB0,marginal.likelihood="Chib95")
	# quadratic
	m2 <- MCMCregress(y~x+I(x^2),data=mydata,c0=cc0,d0=dd0,B0=BB0,marginal.likelihood="Chib95")
	BF <- BayesFactor(m0,m1,m2)
	PPM <- PostProbMod(BF)

	list(m0=m0,m1=m1,m2=m2,BF=BF,PPM=PPM)
}

A <- Bregs()
A$P

####################
#  plot a sample of posterior curves
ppc <- function(mcmcout,x,k=100,type="p1") {
	ri <- sample(1:dim(mcmcout)[1],k)
	S <- mcmcout[ri,]
	if (type == "p1") { # y = a + bx
		for (j in 1:k) abline(S[j,1:2],col="pink")
		abline(mean(mcmcout[,1]),mean(mcmcout[,2]),col="red")
		points(x,attr(mcmcout,"y"),col="red")
	}
	if (type == "p0") { # y = a
		for (j in 1:k) abline(h=S[j,1],col="blue",add=TRUE)
		a <- c(mean(mcmcout[,1]))
		abline(h=a,col="black",add=TRUE)
	}
	if (type =="p2") { # y = a0+a1x+a2x^2
		for (j in 1:k) curve(S[j,1]+S[j,2]*x+S[j,3]*x^2,col="brown",add=TRUE)
		b <- colMeans(mcmcout)
		curve(b[1]+b[2]*x+b[3]*x^2,col="black",add=TRUE)
	}
}

plot(x,y)
ppc(m0,x,type="p0")

plot(x,y)
ppc(m1,x,type="p1")

plot(x,y)
ppc(m2,x,type="p2")

#####################################################################
# non-linear regression.
# read a data.frame from the web:
sample1 <- read.table("http://acccn.net/cr569/Rstuff/mat308/yearpop.txt",header=T)
abline(lm(Population~Year,data=sample1)$coef,col="red")
p1 <- lm(Population~Year,data=sample1)
plot(p1)
y <- sample1$Pop
x <- sample1$Year - 1960
plot(x,y)

# quadratic poly:
p2 <- lm(y~x+I(x^2))
a <- p2$coef
curve(a[1]+a[2]*x+a[3]*x^2,col="red",add=T)
plot(p2)

# cubic:
plot(x,y)
p3 <- lm(y~x+I(x^2)+I(x^3))
a <- p3$coef
curve(a[1]+a[2]*x+a[3]*x^2+a[4]*x^3,add=T,col="blue")

# kpoly:
kpoly <- function(k=2,x,y,...) {
    f <- "1"
    for (j in 1:k) f <- paste(f,"+I(x^",j,")",sep="")
    pk <- lm(as.formula(paste("y~",f,sep="")))
    a <- pk$coef
    h <- "a[1]"
    for (j in 1:k) h <- paste(h,"+a[",j,"]*x^",j,sep="")
    plot(x,y)
    curve(h,add=T,...)
    pk
}