## Bayes h2018: Exam project 2012, Exercise 2 dd <- read.table("cars2012.txt",sep=" ",header=T) #the data from the exercise n <- dim(dd)[1] ### a) plot(dd[,"x"],dd[,"y"],xlab="years",ylab="precent of new price",ylim=c(0,100),xlim=c(1,10)) olm <- lm(y ~ x, data=dd) lines(dd[,"x"],olm$fitted.values,type="l",lwd=3,col="red") ### b) Study reasonable prior values: sim <- 10^3 asim <- runif(sim,80,100) bsim <- rgamma(sim,0.5,1) hist(bsim) csim <- rgamma(sim,20,1) hist(csim) at <- 95 bt <- 0.3 ct <- 10 epssim <- rgamma(n,ct,ct) ysim <- at*exp(-bt*dd[,"x"])*epssim plot(dd[,"x"],ysim,xlab="years",ylab="precent of new price",ylim=c(0,100),xlim=c(1,10)) ### c) ML mlln <- function(param,yy,xx){ a <- param[1] b <- param[2] c <- param[3] -sum(log(dgamma(yy/(a*exp(-b*xx)),c,c)) - log(a*exp(-b*xx))) } cc <- optim(c(85,0.3,44),mlln,yy=dd[,"y"],xx=dd[,"x"]) xval <- seq(1,10,length=100) plot(dd[,"x"],dd[,"y"],xlab="years",ylab="precent of new price",ylim=c(0,100),xlim=c(1,10)) olm <- lm(y ~ x, data=dd) lines(dd[,"x"],olm$fitted.values,type="l",lwd=3,col="red") lines(xval,cc$par[1]*exp(-cc$par[2]*xval),col="blue",lwd=3) ### d) Lazy Bayes: cc <- optim(c(90,0.3,44),mlln,yy=dd[,"y"],xx=dd[,"x"],hessian=T) covM <- solve(cc$hessian) round(cov2cor(covM),4) round(cc$par,3) # 90.897; 0.285; 44.071 round(c(cc$par[1] - 1.645*sqrt(covM[1,1]),cc$par[1] + 1.645*sqrt(covM[1,1])),3) #(80.076;101.717) round(c(cc$par[2] - 1.645*sqrt(covM[2,2]),cc$par[2] + 1.645*sqrt(covM[2,2])),3) #(0.263;0.307) round(c(cc$par[3] - 1.645*sqrt(covM[3,3]),cc$par[3] + 1.645*sqrt(covM[3,3])),3) #(25.423;62.719) ## MCMC lupost <- function(param,yy,xx){ a <- param[1] b <- param[2] c <- param[3] if( c < 0) return(-Inf) -mlln(c(a,b,c),yy,xx) + log(dunif(a,30,100)) + log(dunif(b,0.1,5)) + log(dunif(c,0.5,90)) } sim <- 10^4 dsim <- matrix(NA,sim,3) aprob <- rep(NA,sim) dsim[1,] <- cc$par # starting in MLE for (i in 2:sim){ dprop <- dsim[i-1,] + 1*sqrt(diag(covM))*rnorm(3) aprob[i] <- min(1,exp(lupost(dprop,yy=dd[,"y"],dd[,"x"]) - lupost(dsim[i-1,],yy=dd[,"y"],dd[,"x"]))) keep <- runif(1) <= aprob[i] dsim[i,] <- keep*dprop + (1-keep)*dsim[i-1,] } mean(aprob,na.rm=T) dss <- dsim[(sim/2):sim,] plot(dss[,3]) plot(dss[,2]) plot(dss[,1]) colMeans(dss) cor(dss) round(quantile(dss[,1],c(0.05,0.5,0.95)),3) #(80.828;98.894 ) round(quantile(dss[,2],c(0.05,0.5,0.95)),3) #(0.263;0.303) round(quantile(dss[,3],c(0.05,0.5,0.95)),3) #(26.304;63.777) ### e) Half time: t0sim <- (log(dss[,1])-log(50))/dss[,2] hist(t0sim) round(quantile(t0sim,c(0.05,0.5,0.95)),3) #(1.794;2.281) ### f) Predict: # Expected price: expsim <- dss[,1]*exp(-dss[,2]*5) hist(expsim) round(quantile(expsim,c(0.05,0.5,0.95)),3) #(20.882;22.900) # Posterior predictive: predsim <- dss[,1]*exp(-dss[,2]*5)*rgamma(5001,dss[,3],dss[,3]) hist(predsim) round(quantile(predsim,c(0.05,0.5,0.95)),3) #(16.327;27.853) plot(density(expsim,adjust=2),lwd=2,xlim=c(10,35)) lines(density(predsim,adjust=1.5),col="red",lwd=2)