pausa <- function() { cat ("\n >> press Enter to continue\n"); scan() }

m.true = 2; c.true = 1; sigma.true=2 
x = 1:20
y = m.true*x + c.true + rnorm(length(x), 0, sigma.true) 

plot(x,y, col='blue',ylim=c(0,max(y)) )
pausa()

data = list(x=x, y=y, N=length(x))
inits = list(m=1, c=0, tau=1)

model = "tmp_model.bug"
write("
var mu.y[N];
model{
    for (i in 1:N) {
      y[i] ~ dnorm(mu.y[i], tau);
      mu.y[i] <- x[i]*m + c;
    }
    c ~ dnorm(0, 1.0E-6);
    m ~ dnorm(0, 1.0E-6);
    tau ~ dgamma(1.0, 1.0E-6);
    sigma <- 1.0/sqrt(tau); 
}
", model)

library(rjags)

ns=10000
jm <- jags.model(model, data, inits)
update(jm, 100)
chain <- coda.samples(jm, c("c","m","sigma"), n.iter=ns)

plot(chain)
print(summary(chain))
pausa()

c     <- as.vector(chain[[1]][,1])
m     <- as.vector(chain[[1]][,2])
sigma <- as.vector(chain[[1]][,3])
plot(x,y, col='blue',ylim=c(0,max(y)) )
abline(mean(c), mean(m), col='red')
pausa()

abline(lm(y~x), col='black')   # least squares
pausa()

plot(m,c,col='cyan')
cat(sprintf("rho(m,x) = %.3f\n", cor(m,c) ))
pausa()

# estrapolation
xf = 30
yf = m*xf + c + rnorm(ns, 0, sigma)

hist(yf, nc=100, col='magenta')
cat(sprintf("y(x=%.0f) = %.2f +- %.2f\n", xf, mean(yf), sd(yf)) )