------------------------------ FILE: nonlinearmodel-stan.stan ---------------------------
data {
  int<lower=1> n;
  real x[n];
  real y[n];
}

parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}

model {
  alpha ~ normal(2,1);
  beta  ~ normal(3,1);
  sigma ~ inv_gamma(10,0.018);
  for (i in 1:n)
   y[i] ~ normal(alpha*x[i]/(beta+x[i]),sigma);
}



------------------------------ FILE: nonlinearmodel-stan.R ---------------------------

#install.packages("rstan")
library("rstan")
options(mc.cores=4)

# Data
substrate = c(0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 15.0)
rate      = c(0.26, 0.50, 0.85, 1.10, 1.40, 1.55, 1.65, 1.75)
x         = substrate
y         = rate
n         = length(x)

# Running stan code
model = stan_model("nonlinearmodel-stan.stan")

fit = sampling(model,list(n=n,y=y,x=x),iter=35000,warmup=10000,thin=10,chains=4)

print(fit)

params = extract(fit)

par(mfrow=c(2,3))
ts.plot(params$alpha,xlab="Iterations",ylab="alpha")
ts.plot(params$beta,xlab="Iterations",ylab="beta")
ts.plot(params$sigma,xlab="Iterations",ylab="sigma")
hist(params$alpha,main="",xlab="alpha")
hist(params$beta,main="",xlab="beta")
hist(params$sigma,main="",xlab="sigma")

par(mfrow=c(1,3))
plot(params$alpha,params$beta)
plot(params$alpha,params$sigma)
plot(params$beta,params$sigma)


alpha = params$alpha
beta  = params$beta
sigma = params$sigma
M     = length(alpha)
N     = 20
xx    = seq(min(x),max(x),length=N)
yy    = matrix(0,N,M)
for (i in 1:N)
  yy[i,] = rnorm(M,alpha*xx[i]/(beta+xx[i]),sigma)
quant = t(apply(yy,1,quantile,c(0.05,0.5,0.95)))

par(mfrow=c(1,1))
plot(x,y,ylab="Rate of a reaction (μmol)",xlab="Concentration of a substrate (mM)",
     ylim=range(quant,y),pch=16,cex=0.75)
lines(xx,quant[,1],lty=2,col=2,lwd=2)
lines(xx,quant[,2],col=2,lwd=2)
lines(xx,quant[,3],lty=2,col=2,lwd=2)
points(x,y,cex=0.75,pch=16)
N        = 200
yy       = seq(1.5,1.8,length=N)
postpred = rep(0,N)
for (i in 1:N)
  postpred[i] = mean(dnorm(yy[i],alpha*10/(beta+10),sigma))

par(mfrow=c(1,1))
plot(yy,postpred,ylab="Posterir predictive density",xlab="Rate of a reaction (μmol)",
     type="l",col=2,ylim=c(0,15),lwd=2)
legend("topright",legend="x = 10",bty="n")
points(1.65,0,pch=16,cex=2)