if (FALSE) {
#Load the library
library(brnn)
###############################################################
#Example 1
#Noise triangle wave function, similar to example 1 in Foresee and Hagan (1997)
#Generating the data
x1=seq(0,0.23,length.out=25)
y1=4*x1+rnorm(25,sd=0.1)
x2=seq(0.25,0.75,length.out=50)
y2=2-4*x2+rnorm(50,sd=0.1)
x3=seq(0.77,1,length.out=25)
y3=4*x3-4+rnorm(25,sd=0.1)
x=c(x1,x2,x3)
y=c(y1,y2,y3)
#With the formula interface
out=brnn(y~x,neurons=2)
#With the default S3 method the call is
#out=brnn(y=y,x=as.matrix(x),neurons=2)
plot(x,y,xlim=c(0,1),ylim=c(-1.5,1.5),
main="Bayesian Regularization for ANN 1-2-1")
lines(x,predict(out),col="blue",lty=2)
legend("topright",legend="Fitted model",col="blue",lty=2,bty="n")
###############################################################
#Example 2
#sin wave function, example in the Matlab 2010b demo.
x = seq(-1,0.5,length.out=100)
y = sin(2*pi*x)+rnorm(length(x),sd=0.1)
#With the formula interface
out=brnn(y~x,neurons=3)
#With the default method the call is
#out=brnn(y=y,x=as.matrix(x),neurons=3)
plot(x,y)
lines(x,predict(out),col="blue",lty=2)
legend("bottomright",legend="Fitted model",col="blue",lty=2,bty="n")
###############################################################
#Example 3
#2 Inputs and 1 output
#the data used in Paciorek and
#Schervish (2004). The data is from a two input one output function with Gaussian noise
#with mean zero and standard deviation 0.25
data(twoinput)
#Formula interface
out=brnn(y~x1+x2,data=twoinput,neurons=10)
#With the default S3 method
#out=brnn(y=as.vector(twoinput$y),x=as.matrix(cbind(twoinput$x1,twoinput$x2)),neurons=10)
f=function(x1,x2) predict(out,cbind(x1,x2))
x1=seq(min(twoinput$x1),max(twoinput$x1),length.out=50)
x2=seq(min(twoinput$x2),max(twoinput$x2),length.out=50)
z=outer(x1,x2,f) # calculating the density values
transformation_matrix=persp(x1, x2, z,
main="Fitted model",
sub=expression(y==italic(g)~(bold(x))+e),
col="lightgreen",theta=30, phi=20,r=50,
d=0.1,expand=0.5,ltheta=90, lphi=180,
shade=0.75, ticktype="detailed",nticks=5)
points(trans3d(twoinput$x1,twoinput$x2, f(twoinput$x1,twoinput$x2),
transformation_matrix), col = "red")
###############################################################
#Example 4
#Gianola et al. (2011).
#Warning, it will take a while
#Load the Jersey dataset
data(Jersey)
#Fit the model with the FULL DATA
#Formula interface
out=brnn(pheno$yield_devMilk~G,neurons=2,verbose=TRUE)
#Obtain predictions and plot them against fitted values
plot(pheno$yield_devMilk,predict(out))
#Predictive power of the model using the SECOND set for 10 fold CROSS-VALIDATION
data=pheno
data$X=G
data$partitions=partitions
#Fit the model for the TESTING DATA
out=brnn(yield_devMilk~X,
data=subset(data,partitions!=2),neurons=2,verbose=TRUE)
#Plot the results
#Predicted vs observed values for the training set
par(mfrow=c(2,1))
plot(out$y,predict(out),xlab=expression(hat(y)),ylab="y")
cor(out$y,predict(out))
#Predicted vs observed values for the testing set
yhat_R_testing=predict(out,newdata=subset(data,partitions==2))
ytesting=pheno$yield_devMilk[partitions==2]
plot(ytesting,yhat_R_testing,xlab=expression(hat(y)),ylab="y")
cor(ytesting,yhat_R_testing)
}
Run the code above in your browser using DataLab