#################### COMANDOS DE WINBUGS ############################

#EJERCICIOS DE CLASE: 

#Ventanas necesarias: 
#1) Model -> Specification
#2) Model -> Update
#3) Inference -> Samples
#4) Sample Monitor Tool (coda)


######################################################################
#Ejercicio 1

#---------------------------------------------------------------------
#(a)

model
{
for (i in 1:N) {
	credito[i] ~ dbin(theta,1);
	}
theta ~ dbeta(89.4,59.6);
phi <- theta/(1-theta);
}

list(
	credito = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
	N=100)

list(theta=0.5)

#--------------------------------------------------------------------
#(b)

model
{
for (i in 1:N) {
	credito[i] ~ dbin(theta,1);
	}
theta <- exp(lambda)/(1+exp(lambda));
lambda ~ dnorm(0.4054,33.33);
phi <- theta/(1-theta);

}

list(
	credito = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
	N=100)

list(lambda=0.4)

#--------------------------------------------------------------------
#(c)

model
{
for (i in 1:N) {
	credito[i] ~ dbin(theta,1);
	}
theta ~ dbeta(0.5,0.5);
phi <- theta/(1-theta);
}

list(
	credito = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
	N=100)

list(theta=0.5)


######################################################################
#Ejercicio 2

#---------------------------------------------------------------------
#(a)

model
{
for (i in 1:N) {
	utilidad[i] ~ dnorm(mu,tau);
	}
mu ~ dnorm(200,0.025);
tau <- 1/(sig*sig);
sig ~ dgamma(10,1);
utilidadF ~ dnorm(mu,tau);
}

list(
	utilidad = c(212, 207, 210, 196, 223, 193, 196, 210, 202, 221),
	N=10)

list(mu=200, sig=10, utilidadF=200)

#---------------------------------------------------------------------
#(b)

model
{
for (i in 1:N) {
	utilidad[i] ~ dnorm(mu,tau);
	}
mu ~ dnorm(0,0.00001);
tau ~ dgamma(0.0001,0.0001);
utilidadF ~ dnorm(mu,tau);
}

list(
	utilidad = c(212, 207, 210, 196, 223, 193, 196, 210, 202, 221),
	N=10)

list(mu=0, tau=1, utilidadF=0)


######################################################################
#Ejercicio 3

model
{
for (i in 1:N) {
	SP[i] ~ dnorm(mu[i],tau);
	mu[i] <- beta[1]+beta[2]*MO[i]
	}
for (i in 1:2) {
	beta[i] ~ dnorm(0.0,0.0001);
	}
tau ~ dgamma(0.0001,0.0001);
for (i in 1:N) {
	mu1[i] <- beta[1]+beta[2]*MO[i];
	SP1[i] ~ dnorm(mu1[i],tau)
	}
mu2 <- beta[1]+beta[2]*MO2;
SP2 ~ dnorm(mu2,tau);
}


list(
	SP = c(3.1,2.3,3,1.9,2.5,3.7,3.4,2.6,2.8,1.6,2,2.9,2.3,3.2,1.8,1.4,2,3.8,2.2,1.5),
	MO = c(5.5,4.8,4.7,3.9,4.5,6.2,6,5.2,4.7,4.3,4.9,5.4,5,6.3,4.6,4.3,5,5.9,4.1,4.7),
	N=20, MO2=3)

list(beta=c(0,0),tau=1,SP2=0)


######################################################################
#Ejercicio 4

model
{
for (i in 1:N) {
	Y[i] ~ dnorm(mu[i],tau);
	mu[i] <- alpha+beta[1]*X1[i]+beta[2]*X2[i]+beta[3]*X3[i];
	}
for (i in 1:3) {
	beta[i] ~ dnorm(0.0,0.0001);
	}
alpha ~ dnorm(0.0,0.0001);
tau ~ dgamma(0.0001,0.0001);
for (i in 1:N) {
	Y1[i] ~ dnorm(mu1[i],tau);
	mu1[i] <- alpha+beta[1]*X1[i]+beta[2]*X2[i]+beta[3]*X3[i];
	}
for (i in 1:3) {
	Y2F[i] ~ dnorm(mu2[i],tau);
	mu2[i] <- alpha+beta[1]*X1F[i]+beta[2]*X2F[i]+beta[3]*X3F[i];
	}
}


list(
	Y=c(33.2, 40.3, 38.7, 46.8, 41.4, 37.5, 39., 40.7, 30.1, 52.9, 38.2, 31.8, 
	  43.3, 44.1, 42.8, 33.6, 34.2, 48., 38., 35.9, 40.4, 36.8, 45.2, 35.1), 
	X1=c(3.5, 5.3, 5.1, 5.8, 4.2, 6., 6.8, 5.5, 3.1, 7.2, 4.5, 4.9, 8., 6.5, 6.6, 
	  3.7, 6.2, 7., 4., 4.5, 5.9, 5.6, 4.8, 3.9), 
	X2=c(9., 20., 18., 33., 31., 13., 25., 30., 5., 47., 25., 11., 23., 35., 39., 
	  21., 7., 40., 35., 23., 33., 27., 34., 15.), 
	X3=c(6.1, 6.4, 7.4, 6.7, 7.5, 5.9, 6., 4., 5.8, 8.3, 5., 6.4, 7.6, 7., 5., 4.4, 
	  5.5, 7., 6., 3.5, 4.9, 4.3, 8., 5.),
	X1F=c(5.4, 6.2, 6.4), 
	X2F=c(17, 12, 21),
	X3F=c(6.0, 5.8, 6.1),
	N=24)

list(alpha=0,beta=c(0,0,0),tau=1,Y2F=c(0,0,0))


######################################################################
#Ejercicio 5

#---------------------------------------------------------------------
#(a) Modelo estático (x[t])

model  {
	tau.y  <- 1
	for (t in 1:13){y[t] ~ dnorm(mu[t],tau.y);
                          mu[t] <-  alpha + x[t]*beta;}
# forecast
	for (t in 1:13){  mu.new[t] <- alpha + x[t]*beta;
                       y.new[t] ~ dnorm(mu.new[t],tau.y) }
# initials
	alpha ~ dnorm(10,0.001);
	beta ~ dnorm(10,0.001) 
}


list(
	y=c(117., 118.6, 120., 115.5, 115.6, 115.4, 120.2, 122.7, 121.5, 123.4, 128.5, 
	  130., 135.8), 
	x=c(12., 11.8, 11.7, 11.4, 11.2, 11.1, 11., 11., 10.8, 10.7, 10.8, 10.9, 11.)) 


list(alpha=0, beta=10)

#---------------------------------------------------------------------
#(b) Modelo estático (Tiempo)

model  {
	tau.y  <- 1
	for (t in 1:13){y[t] ~ dnorm(mu[t],tau.y);
                          mu[t] <-  alpha + tiempo[t]*beta;}
# forecast
	for (t in 1:13){  mu.new[t] <- alpha + tiempo[t]*beta;
                       y.new[t] ~ dnorm(mu.new[t],tau.y) }
# initials
	alpha ~ dnorm(0,0.001);
	beta ~ dnorm(0,0.001) 
}


list(
	y=c(117., 118.6, 120., 115.5, 115.6, 115.4, 120.2, 122.7, 121.5, 123.4, 128.5, 
	  130., 135.8), 
	tiempo=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)) 


list(alpha=0, beta=0)

#---------------------------------------------------------------------
#(c) Modelo estático (Tiempo + x[t])

model  {
	tau.y  <- 1
	for (t in 1:13){y[t] ~ dnorm(mu[t],tau.y);
                          mu[t] <-  alpha + beta[1]*tiempo[t] + beta[2]*x[t];}
# forecast
	for (t in 1:13){  mu.new[t] <- alpha + beta[1]*tiempo[t] + beta[2]*x[t];
                       y.new[t] ~ dnorm(mu.new[t],tau.y) }
# initials
	alpha ~ dnorm(0,0.001);
	beta[1] ~ dnorm(0,0.001);
	beta[2] ~ dnorm(0,0.001)
}


list(
	y=c(117., 118.6, 120., 115.5, 115.6, 115.4, 120.2, 122.7, 121.5, 123.4, 128.5, 
	  130., 135.8), 
	tiempo=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13),
	x=c(12., 11.8, 11.7, 11.4, 11.2, 11.1, 11., 11., 10.8, 10.7, 10.8, 10.9, 11.)) 


list(alpha=0, beta=c(0,0))

#---------------------------------------------------------------------
#(d) Modelo dinámico (con variable explicativa)

model  {# fixed precisions, as in Harrison and West, Ch 3
   W <- 0.05; W.inv <- 1/W
   tau.y  <- 1
#  W.inv ~ dgamma(0.01,0.01);
#  tau.y ~ dgamma(0.01,0.01);

# observation model 
     for (t in 1:13){y[t] ~ dnorm(mu[t],tau.y);
                          mu[t] <-  x[t]*beta[t];}
# state model
     for (t in 2:13){ beta[t]  ~ dnorm(beta[t-1],W.inv);}
# one-step ahead forecasts
     for (t in 2:13){  beta.new[t] ~ dnorm(beta[t-1],W.inv);
		       mu.new[t] <- beta.new[t]*x[t];
                       y.new[t] ~ dnorm(mu.new[t],tau.y) }
# settings for year 1
      beta[1] ~ dnorm(10,0.001) 
}


list(
	y=c(117., 118.6, 120., 115.5, 115.6, 115.4, 120.2, 122.7, 121.5, 123.4, 128.5, 
	  130., 135.8), 
	x=c(12., 11.8, 11.7, 11.4, 11.2, 11.1, 11., 11., 10.8, 10.7, 10.8, 10.9, 11.)) 


list(beta=c(10,10,10,10,10,10,10,10,10,10,10,10,10))

#---------------------------------------------------------------------
#(e) Modelo dinámico (sin variables explicativas)

model  {# fixed precisions, as in Harrison and West, Ch 3
#   W <- 10; W.inv <- 1/W
#   tau.y  <- 0.01 
    W.inv ~ dgamma(0.01,0.01);
    tau.y ~ dgamma(0.01,0.01);

# observation model 
     for (t in 1:13){y[t] ~ dnorm(mu[t],tau.y);
		     mu[t] <- beta[t]}
# state model
     for (t in 2:13){ beta[t] ~ dnorm(beta[t-1],W.inv);}
# one-step ahead forecasts
     for (t in 2:13){  beta.new[t] ~ dnorm(beta[t-1],W.inv);
		       mu.new[t] <- beta.new[t];
                       y.new[t] ~ dnorm(mu.new[t],tau.y);}
# settings for year 1
      beta[1] ~ dnorm(0,0.001) 
}


list(
	y=c(117., 118.6, 120., 115.5, 115.6, 115.4, 120.2, 122.7, 121.5, 123.4, 128.5, 
	  130., 135.8))

list(beta=c(0,0,0,0,0,0,0,0,0,0,0,0,0))


######################################################################
#Ejercicio 6

#---------------------------------------------------------------------
#(a) Modelo estático

model {for (i in 1:N) {SHARE[i] ~ dnorm(mu[i],P)
mu[i] <- b[1]+b[2]*PRICE[i]+b[3]*OPROM[i]+b[4]*CPROM[i]}
for (i in 2:N) {SHARE.STEP[i] ~ dnorm(mu[i-1],P)
                     resid.step[i] <- SHARE.STEP[i] - SHARE[i]
                     absdev[i] <- abs(resid.step[i])}
MAD <- sum(absdev[2:N])/(N-1)
P ~ dgamma(0.5,0.5)
b[1] ~ dnorm(42,0.04);    
for(j in 2:4){b[j] ~ dnorm(0,0.25)}}


list(N=104,
CPROM=c(0.74, 0.758, 0.769, 0.781, 0.833, 0.841, 0.806, 0.825, 0.849, 0.859, 0.816, 0.67, 
  0.575, 0.529, 0.453, 0.385, 0.326, 0.275, 0.148, 0.032, -0.145, -0.302, -0.46, -0.609, 
  -0.738, -0.817, -0.849, -0.846, -0.837, -0.846, -0.862,	-0.85, -0.807, -0.749, -0.768, 
  -0.853, -0.923, -1.003, -1.1, -1.183, -1.22, -1.158, -1.008, -0.818, -0.561, -0.263, 
  0.031, 0.326, 0.606, 0.846, 1.018, 1.079, 1.013, 0.827, 0.539, 0.215, -0.086, -0.346, 
  -0.543, -0.641,	-0.646, -0.597, -0.541, -0.514, -0.528, -0.591, -0.709, -0.841, -0.935, 
  -1.01, -1.065, -1.085, -1.093, -1.132, -1.224, -1.348, -1.459, -1.539, -1.555, -1.504, 
  -1.37, -1.141, -0.857, -0.514, -0.134, 0.216, 0.53, 0.842, 1.145, 1.405, 1.623, 1.809, 
  1.932, 1.998, 2.033, 2.036, 1.986, 1.851, 1.657, 1.478, 1.35, 1.24, 1.119, 1.056), 
PRICE=c(-0.463, -0.607, -0.724, -0.83, -0.937, -0.961, -0.938, -0.883, -0.748, -0.574, 
  -0.351, -0.106, 0.11, 0.337, 0.562, 0.744, 0.945, 1.153, 1.383, 1.56, 1.689, 1.77, 
  1.881, 1.907, 2.008, 2.036, 2.133, 2.212, 2.289, 2.292, 2.26, 2.107, 1.945, 1.818, 
  1.644, 1.479, 1.322, 1.215, 1.105, 1.024, 0.893, 0.802, 0.717, 0.705, 0.648, 0.636, 
  0.598, 0.609, 0.585, 0.534, 0.47, 0.412, 0.312, 0.167, -0.03, -0.223, -0.365, -0.516, 
  -0.674, -0.795, -0.866, -0.91, -0.871, -0.787, -0.681, -0.531, -0.49, -0.508, -0.548, 
  -0.668, -0.856, -1.055, -1.23, -1.452, -1.588, -1.662, -1.655, -1.61, -1.558, -1.481, 
  -1.439, -1.348, -1.307, -1.311, -1.33, -1.391, -1.372, -1.329, -1.317, -1.28, -1.228, 
  -1.178, -1.106, -1.071, -1.016, -0.881, -0.736, -0.615, -0.471, -0.265, -0.083, 0.052, 
  0.132, 0.206), 
OPROM=c(-0.11, -0.269, -0.425, -0.487, -0.596, -0.456, -0.161, 0.117, 0.537, 0.934, 
  1.297, 1.588, 1.79, 1.917, 1.992, 1.953, 1.822, 1.566, 1.287, 0.904, 0.517, 0.273, 
  -0.029, -0.12, 0.064, 0.281, 0.521, 0.75, 0.826, 0.864, 0.923, 0.802, 0.604, 0.396, 
  0.161, -0.134, -0.367, -0.651, -0.841, -0.976, -1.044, -1.154, -1.174, -1.192, -1.216, 
  -1.201, -1.171, -1.153, -1.153, -1.154, -1.204, -1.201, -1.196, -1.211, -1.258, -1.271, 
  -1.2, -1.048, -0.83, -0.467, -0.108, 0.341, 0.746, 1.047, 1.342, 1.618, 1.652, 1.613, 
  1.325, 1.079, 0.769, 0.399, -0.025, -0.367, -0.665, -0.777, -0.957, -0.994, -1.049, 
  -1.063, -1.113, -1.153, -1.178, -1.209, -1.209, -1.208, -1.241, -1.272, -1.273, -1.297, 
  -1.292, -1.279, -1.128, -0.836, -0.478, 0.028, 0.545, 1.14, 1.59, 1.975, 2.2, 2.387, 
  2.442, 2.423), 
SHARE=c(42.4, 42.4, 42.7, 42.8, 42.3, 42.5, 42.7, 42.1, 42.1, 41.4, 41.4, 41.8, 41.3, 
  41.3, 41.2, 41.1, 40.8, 40.3, 39.8, 40.2, 39.6, 39.8, 39.5, 39.8, 39.8, 39.8, 39.7, 
  39.4, 39.4, 39.4, 39.7, 39.7, 39.8, 39.2, 40.1, 40.1, 40.2, 40.3, 40.7, 40.8, 41.3, 
  40.9, 41., 41.1, 40.7, 40.6, 40.6, 40.4, 40.3, 40.1, 39.9, 40., 40.2, 40.5, 41.4, 41.7, 
  42.1, 42.4, 42.9, 42.9, 43.1, 42.8, 43., 42.7, 42.8, 42.5, 43., 43., 42.9, 43., 43.2, 
  43.1, 43.2, 43.3, 43.6, 43.8, 43.8, 44., 43.6, 43.9, 43.5, 43.3, 43.3, 43., 42.7, 43., 
  42.9, 42.7, 42.5, 42.6, 42.5, 42.3, 42.4, 42.5, 42.4, 42.3, 42.2, 42.5, 42.7, 42.5, 
  42.1, 42.1, 42.3, 41.9))


list(b=c(42,0,0,0),P=1)


#---------------------------------------------------------------------
#(b) Modelo dinámico

model   {     # observational variance with discounting
                    P ~ dgamma(0.5,0.5)
                    Pobs[1] <- P
   for (t in 2:N) { Pobs[t] <- P*pow(0.99,t-1)}
# one step ahead forecasts
   for (t in 2:N) { SHARE.one[t]   ~ dnorm(mu[t-1],Pobs[t-1])
                        absdev[t] <- abs(SHARE.one[t] - SHARE[t])}
                        MAD <- sum(absdev[2:N])/(N-1)
# regression priors with discounting
      for (j in 1:4) {Pbeta[j] ~ dgamma(0.1,0.1)}
      for (t in 2:N) {b[t,1] ~ dnorm(b[t-1,1],Pb[t,1])
                      b[t,2] ~ dnorm(b[t-1,2],Pb[t,2]) I(,0)
                      b[t,3] ~ dnorm(b[t-1,3],Pb[t,3]) I(0,)
                      b[t,4] ~ dnorm(b[t-1,4],Pb[t,4]) I(,0)
# discount factors for level and predictors
                      DF[1,t] <- pow(0.99,t-2)
                      DF[2,t] <- pow(1,t-2)
                      DF[3,t] <- pow(1,t-2)
                      DF[4,t] <- pow(1,t-2)
      for (j in 1:4) {Pb[t,j] <- DF[j,t]*Pbeta[j]}   
}
# Initialising regression priors
      b[1,1] ~ dnorm(42,0.04)
      b[1,2] ~ dnorm(0,0.25)  I(,0)
      b[1,3] ~ dnorm(0,0.25)  I(0,)
      b[1,4] ~ dnorm(0,0.25)  I(,0)
# on-line mean 
      for (t in 1:N) {  SHARE[t]   ~ dnorm(mu[t],Pobs[t])
               mu[t] <-    b[t,1] + b[t,2]*PRICE[t]  + b[t,3]*OPROM[t] + b[t,4]*CPROM[t]} }


list(N=104,
CPROM=c(0.74, 0.758, 0.769, 0.781, 0.833, 0.841, 0.806, 0.825, 0.849, 0.859, 0.816, 0.67, 
  0.575, 0.529, 0.453, 0.385, 0.326, 0.275, 0.148, 0.032, -0.145, -0.302, -0.46, -0.609, 
  -0.738, -0.817, -0.849, -0.846, -0.837, -0.846, -0.862,	-0.85, -0.807, -0.749, -0.768, 
  -0.853, -0.923, -1.003, -1.1, -1.183, -1.22, -1.158, -1.008, -0.818, -0.561, -0.263, 
  0.031, 0.326, 0.606, 0.846, 1.018, 1.079, 1.013, 0.827, 0.539, 0.215, -0.086, -0.346, 
  -0.543, -0.641,	-0.646, -0.597, -0.541, -0.514, -0.528, -0.591, -0.709, -0.841, -0.935, 
  -1.01, -1.065, -1.085, -1.093, -1.132, -1.224, -1.348, -1.459, -1.539, -1.555, -1.504, 
  -1.37, -1.141, -0.857, -0.514, -0.134, 0.216, 0.53, 0.842, 1.145, 1.405, 1.623, 1.809, 
  1.932, 1.998, 2.033, 2.036, 1.986, 1.851, 1.657, 1.478, 1.35, 1.24, 1.119, 1.056), 
PRICE=c(-0.463, -0.607, -0.724, -0.83, -0.937, -0.961, -0.938, -0.883, -0.748, -0.574, 
  -0.351, -0.106, 0.11, 0.337, 0.562, 0.744, 0.945, 1.153, 1.383, 1.56, 1.689, 1.77, 
  1.881, 1.907, 2.008, 2.036, 2.133, 2.212, 2.289, 2.292, 2.26, 2.107, 1.945, 1.818, 
  1.644, 1.479, 1.322, 1.215, 1.105, 1.024, 0.893, 0.802, 0.717, 0.705, 0.648, 0.636, 
  0.598, 0.609, 0.585, 0.534, 0.47, 0.412, 0.312, 0.167, -0.03, -0.223, -0.365, -0.516, 
  -0.674, -0.795, -0.866, -0.91, -0.871, -0.787, -0.681, -0.531, -0.49, -0.508, -0.548, 
  -0.668, -0.856, -1.055, -1.23, -1.452, -1.588, -1.662, -1.655, -1.61, -1.558, -1.481, 
  -1.439, -1.348, -1.307, -1.311, -1.33, -1.391, -1.372, -1.329, -1.317, -1.28, -1.228, 
  -1.178, -1.106, -1.071, -1.016, -0.881, -0.736, -0.615, -0.471, -0.265, -0.083, 0.052, 
  0.132, 0.206), 
OPROM=c(-0.11, -0.269, -0.425, -0.487, -0.596, -0.456, -0.161, 0.117, 0.537, 0.934, 
  1.297, 1.588, 1.79, 1.917, 1.992, 1.953, 1.822, 1.566, 1.287, 0.904, 0.517, 0.273, 
  -0.029, -0.12, 0.064, 0.281, 0.521, 0.75, 0.826, 0.864, 0.923, 0.802, 0.604, 0.396, 
  0.161, -0.134, -0.367, -0.651, -0.841, -0.976, -1.044, -1.154, -1.174, -1.192, -1.216, 
  -1.201, -1.171, -1.153, -1.153, -1.154, -1.204, -1.201, -1.196, -1.211, -1.258, -1.271, 
  -1.2, -1.048, -0.83, -0.467, -0.108, 0.341, 0.746, 1.047, 1.342, 1.618, 1.652, 1.613, 
  1.325, 1.079, 0.769, 0.399, -0.025, -0.367, -0.665, -0.777, -0.957, -0.994, -1.049, 
  -1.063, -1.113, -1.153, -1.178, -1.209, -1.209, -1.208, -1.241, -1.272, -1.273, -1.297, 
  -1.292, -1.279, -1.128, -0.836, -0.478, 0.028, 0.545, 1.14, 1.59, 1.975, 2.2, 2.387, 
  2.442, 2.423), 
SHARE=c(42.4, 42.4, 42.7, 42.8, 42.3, 42.5, 42.7, 42.1, 42.1, 41.4, 41.4, 41.8, 41.3, 
  41.3, 41.2, 41.1, 40.8, 40.3, 39.8, 40.2, 39.6, 39.8, 39.5, 39.8, 39.8, 39.8, 39.7, 
  39.4, 39.4, 39.4, 39.7, 39.7, 39.8, 39.2, 40.1, 40.1, 40.2, 40.3, 40.7, 40.8, 41.3, 
  40.9, 41., 41.1, 40.7, 40.6, 40.6, 40.4, 40.3, 40.1, 39.9, 40., 40.2, 40.5, 41.4, 41.7, 
  42.1, 42.4, 42.9, 42.9, 43.1, 42.8, 43., 42.7, 42.8, 42.5, 43., 43., 42.9, 43., 43.2, 
  43.1, 43.2, 43.3, 43.6, 43.8, 43.8, 44., 43.6, 43.9, 43.5, 43.3, 43.3, 43., 42.7, 43., 
  42.9, 42.7, 42.5, 42.6, 42.5, 42.3, 42.4, 42.5, 42.4, 42.3, 42.2, 42.5, 42.7, 42.5, 
  42.1, 42.1, 42.3, 41.9))


###Usar preferentemente (gen inits)###
list(P=1,Pbeta=c(1,1,1,1),
b=structure(.Data=c(
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,
42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0,42,0,0,0),
.Dim=c(104,4)))


######################################################################
#Ejercicio 7

#---------------------------------------------------------------------
# Modelo Binomial - Varias ligas

model {  
for (i in 1:6) {
	y[i] ~ dbin(pi[i],n[i]);
# alternative form of predictor
#   	z[i] <- log(x[i]+1);
#     	z[i] <- x[i]*x[i];
# alternative links
   	logit(pi[i]) <- beta[1]+beta[2]*x[i];
#   	probit(pi[i]) <- beta[1]+beta[2]*x[i];
#   	cloglog(pi[i]) <- beta[1]+beta[2]*z[i];
# predicted deaths (mean)
	yhat[i] <- n[i]*pi[i];
# goodness of fit (deviance)
   	dev[i] <- y[i]*log(y[i]/yhat[i])+(n[i]-y[i])*log((n[i]-y[i])/(n[i]-yhat[i]))}
  	beta[1] ~ dnorm(0.0,1.0E-3);
  	beta[2] ~ dnorm(0.0,1.0E-3);
# total deviance
  	dev.tot <- 2*sum(dev[]) 
# prediction for a new individual
	y.new ~ dbin(pi.new,n.new);
	logit(pi.new) <- beta[1]+beta[2]*x.new;
#	probit(pi.new) <- beta[1]+beta[2]*x.new;
}


list(
	x=c(0., 5., 30., 75., 150., 250.), 
	y=c(13., 5., 5., 3., 4., 18.), 
	n=c(391., 205., 156., 50., 35., 51.), 
	x.new=200, n.new=100)


list(beta=c(0,0))


######################################################################
#Ejercicio 8

#---------------------------------------------------------------------
# (a) Modelo con tasa en función del tiempo

model {    
for( year in 1 : N ) {
       D[year] ~ dpois(mu[year])    
       log(mu[year]) <- b[1] + b[2]*year }
for (j in 1:2) {b[j] ~ dnorm( 0.0,1.0E-4)}
}


list(N=112,D=c(4,5,4,1,0,4,3,4,0,6,
	3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,
	1,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,
	0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,
	2,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,
	1,0,0,0,0,0,1,0,0,1,0,0))


list(b=c(0,0))


#---------------------------------------------------------------------
# (b) Modelo con tasa constante en dos períodos

model {    
for( year in 1 : N ) {
       D[year] ~ dpois(mu[year])    
       log(mu[year]) <- b[1] + step(year - changeyear) * b[2]    }
for (j in 1:2) {b[j] ~ dnorm( 0.0,1.0E-6)}
actual <- changeyear+1850
changeyear ~ dunif(1,N) 
}


list(N=112,D=c(4,5,4,1,0,4,3,4,0,6,
	3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,
	1,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,
	0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,
	2,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,
	1,0,0,0,0,0,1,0,0,1,0,0))


list(changeyear=50,b=c(0,0))


#---------------------------------------------------------------------
# (c) Modelo dinámico 

model {

for( year in 1 : N ) {
       D[year] ~ dpois(mu[year])
       log(mu[year]) <- b[year] }

# state model
  for (year in 2:N){ b[year] ~ dnorm(b[year-1],W.inv);}

# priors
  b[1] ~ dnorm( 0.0,1.0E-4)
  W.inv ~ dgamma(0.01,0.01)

}


list(N=112,D=c(4,5,4,1,0,4,3,4,0,6,
	3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,
	1,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,
	0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,
	2,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,
	1,0,0,0,0,0,1,0,0,1,0,0))


list(b=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  W.inv=1)


######################################################################
#Ejercicio 9

#---------------------------------------------------------------------
# (a) Modelo con efectos independientes

model { 
for (i in 1:N) {
	incid.rate[Cancer[i],Place[i],Age[i]] ~ dgamma(0.001, 0.001);
	lambda[i] <- incid.rate[Cancer[i],Place[i],Age[i]] * Pops[i]/100000;
	O[i] ~ dpois(lambda[i])	}
}


list(
	O=c(38,13,51,37,5,8,13,20),
	Place=c(1,1,2,2,1,1,2,2),
	Age=c(1,2,1,2,1,2,1,2),
	Cancer=c(1,1,1,1,2,2,2,2),
	Pops=c(103857,155786,135943,203914,103857,155786,135943,203914),
	N=8)


list(incid.rate=structure(.Data=c(20,20,20,20,20,20,20,20),.Dim=c(2,2,2))			


#---------------------------------------------------------------------
# (b) Modelo con efectos intercambiables

model { 
for (i in 1:N) {	
	R[Cancer[i],Place[i],Age[i]] ~ dgamma(alpha, beta);
# rate per million outcome
	lambda1[i] <- R[Cancer[i],Place[i],Age[i]] * Pops[i]/100000
# smr outcome			
	lambda2[i] <- R[Cancer[i],Place[i],Age[i]] * E[i]
# choose outcome
        lambda[i] <- lambda1[i]			
	O[i] ~ dpois(lambda[i])	}
	alpha ~ dexp(1)
	var <- alpha/(beta*beta)
	beta ~ dgamma(1,1)
}

list(
	O=c(38,13,51,37,5,8,13,20),
	E=c(24.08,36.12,31.52,47.28,7.97,11.95,10.43,15.65),
	Place=c(1,1,2,2,1,1,2,2),
	Age=c(1,2,1,2,1,2,1,2),
	Cancer=c(1,1,1,1,2,2,2,2),
	Pops=c(103857,155786,135943,203914,103857,155786,135943,203914),
	N=8)

list(alpha=1,beta=1)



######################################################################
#Ejercicio 10

#---------------------------------------------------------------------
# (a) Modelo con efectos independientes

model 
{
for (j in 1:N) {
	r[j] ~ dbin(pi[j],n[j]);
    	pi[j] ~ dbeta(1,1);}
}


list(n=c(51,16,82,13,43,75,13,10,6,37),
       r=c(24,7,46,9,23,53,8,3,1,23),N=10)


list(pi=c(1,1,1,1,1,1,1,1,1,1))


#---------------------------------------------------------------------
# (b) Modelo con efectos intercambiables

model 
{
for (j in 1:N) {
	r[j] ~ dbin(pi[j],n[j]);
    	pi[j] ~ dbeta(alpha,beta);}
alpha ~ dgamma(1,1);
beta ~ dgamma(1,1);
}


list(n=c(51,16,82,13,43,75,13,10,6,37),
       r=c(24,7,46,9,23,53,8,3,1,23),N=10)


list(alpha=1,beta=1)