## Ajuste de RNA/MLP com R
## para dados simulados: seno
## 22/06/2016

 rm(list=ls(all=TRUE))

 ##- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 ## Exemplo de aplicacao
 ## Cria os dados do exemplo
 set.seed(18)
 X <- seq(0, 2*pi, length.out=90)
 Y <- 3.5*sin(X) + rnorm(length(X), sd=0.8)
 plot(Y ~ X, pch=19, col="blue"); grid()

 library(rjags)

 model_string <- "model{
       # Likelihood
       for(i in 1:N){
           Y[i] ~ dnorm(mu[i], inv.var)
           for(j in 1:H){
               nodo[i,j] <- 2/( 1 + exp( - 2*(w1[j]*x[i] + b1[j]) ) ) - 1
               aux[i,j]  <- nodo[i,j]*w2[j]
           }
           mu[i] <- sum(aux[i,]) + b2
       }

        # Prior for weights
		for( i in 1 : H ) {
			w1[i] ~ dnorm(0, 0.1)
			b1[i] ~ dnorm(0, 0.1)
			w2[i] ~ dnorm(0, 0.1)
		} 
		b2 ~ dnorm(0, 0.01)

       # Prior for inverse variance
		inv.var ~ dgamma(0.01, 0.01)
		sigma <- 1/sqrt(inv.var)
	}"

 # Running the model
 model <- jags.model(textConnection(model_string), data = list(Y=Y, x=X, N=length(Y), H=7), 
                     n.chains = 1, n.adapt= 10000)

 update(model, 10000); # Burnin for 10000 samples

 mcmc_samples <- coda.samples(model, variable.names=c("mu", "sigma"), n.iter=20000)
 summary(mcmc_samples)

 output <- as.matrix(mcmc_samples)

 plot(Y ~ X, pch=19, col="blue"); grid()
 fit    <- apply(output, 2, mean)
 lines(fit[1:length(Y)] ~ X, col="red", lwd=2)

 fit    <- apply(output, 2, median)
 lines(fit[1:length(Y)] ~ X, col="black", lwd=2, lty=2)

 lwr    <- apply(output, 2, function(x) quantile(x, probs=0.025))
 lines(lwr[1:length(Y)] ~ X, col="red", lwd=2, lty=2)
 upr    <- apply(output, 2, function(x) quantile(x, probs=0.975))
 lines(upr[1:length(Y)] ~ X, col="red", lwd=2, lty=2)