library(BARS)
library(gss)
x <- c(1:110)
mu <- 20*abs(exp(-x/100)*cos(x/5))
## replicate example of "bars"
set.seed(20061001,"Mersenne-Twister")
y <- rnorm(110,mu,1)
system.time(out <- bars(x,y,prior="uniform",priorparam=c(1,15),fits=TRUE))
system.time(gssfit1 <- ssanova0(y~x))
pred.gss <- predict(gssfit1,newdata=data.frame(x=x),se=TRUE)
upper <- apply(out$sampfits,2,quantile,0.975)
lower <- apply(out$sampfits,2,quantile,0.025)
upr <- pred.gss$fit+1.96*pred.gss$se.fit
lwr <- pred.gss$fit-1.96*pred.gss$se.fit
bars.out <- sum(mu < lower | upper < mu)
gss.out <- sum(mu < lwr | upr < mu)
layout(matrix(1:2,nrow=1))
m <- min(c(lower,lwr))
M <- max(c(upper,upr))
plot(x,out$postmeans,
type="n",xlab="x",ylab="y",
main=paste("BARS result (",bars.out,"out of",length(x),")"),
ylim=c(m,M))
polygon(c(x,rev(x)),c(lower,rev(upper)),col="grey30",lty=0)
lines(x,mu,col=2,lwd=2)
plot(x,pred.gss$fit,type="n",
xlab="x",ylab="y",
main=paste("GSS result (",gss.out,"out of",length(x),")"),
ylim=c(m,M))
polygon(c(x,rev(x)),c(lwr,rev(upr)),col="grey30",lty=0)
lines(x,mu,col=2,lwd=2)
## replicate example of "barsP"
set.seed(20061001,"Mersenne-Twister")
y <- rpois(110,mu)
system.time(out <- barsP(x,y,prior="uniform",priorparam=c(1,15),fits=TRUE))
system.time(gssfit1 <- gssanova0(y~x,family="poisson"))
pred.gss <- predict(gssfit1,newdata=data.frame(x=x),se=TRUE)
upper <- apply(out$sampfits,2,quantile,0.975)
lower <- apply(out$sampfits,2,quantile,0.025)
upr <- exp(pred.gss$fit+1.96*pred.gss$se.fit)
lwr <- exp(pred.gss$fit-1.96*pred.gss$se.fit)
bars.out <- sum(mu < lower | upper < mu)
gss.out <- sum(mu < lwr | upr < mu)
X11();layout(matrix(1:2,nrow=1))
m <- min(c(lower,lwr))
M <- max(c(upper,upr))
plot(x,out$postmeans,
type="n",xlab="x",ylab="y",
main=paste("BARS result (",bars.out,"out of",length(x),")"),
ylim=c(m,M))
polygon(c(x,rev(x)),c(lower,rev(upper)),col="grey30",lty=0)
lines(x,mu,col=2,lwd=2)
plot(x,pred.gss$fit,type="n",
xlab="x",ylab="y",
main=paste("GSS result (",gss.out,"out of",length(x),")"),
ylim=c(m,M))
polygon(c(x,rev(x)),c(lwr,rev(upr)),col="grey30",lty=0)
lines(x,mu,col=2,lwd=2)