## Start with simple linear model model
library(INLA)

# Generate data

# fixed effect
x = sort(runif(100))
# data vector
y = 1 + 2*x + rnorm(n = 100, sd = 0.1)

# Run inla
formula = y ~  x
result = inla(formula,
              data = list(x = x, y = y),
              family = "gaussian")

# Get summary
summary(result)

# Plot martinal of the intercept
m = result$marginals.fixed[[1]]
plot(m)

# Plot a more smooth function
plot(inla.smarginal(m))

# Extract quantiles
inla.qmarginal(0.05, m)

# Distribution function
inla.pmarginal(0.975, m)

# Density function
inla.dmarginal(1, m)

# Generate 4 realizations
inla.rmarginal(4, m)

# Calculate expected value of x and x^2
E = inla.emarginal(function(x) c(x,x^2), m)
E

# Calculate sd
sqrt(E[2]-E[1]^2)

# Compare to estimate provided by INLA
round(result$summary.fixed[,1:2], 3)

# get the marginal for the precision
tau0 = result$marginals.hyperpar$"Precision for the Gaussian observations"

# Calculate expected value of 1/x and 1/x^2
E = inla.emarginal(function(x) c(1/x,(1/x)^2), tau0)

# Calculate sd
mysd = sqrt(E[2] - E[1]^2)

print(c(mean=E[1], sd=mysd))


## Addition of a random effect

# Generate AR (1) sequence
n = 100
t = 1:n
ar = rep (0 ,n)
for( i in 2:n)
  ar[i] = 0.8 * ar[i -1]+ rnorm( n = 1 , sd = 0.1)
# Generate data with AR (1) component
x = runif(n)
y = 1 + 2 * x + ar + rnorm( n = n , sd = 0.1)
# Run inla
formula = y ~ 1 + x + f(t , model = "ar1" )
result = inla(formula ,
data = list( x = x , y = y , t = t ) ,
    family = "gaussian")
    
summary(result)

## Prediction

# Add new location
x = c(x, 0.3)
t = c(t, 101)
y = c(y, NA)

# Re-compute
result.pred = inla(formula,
                   data = list(x = x, t = t, y = y),
                   family="gaussian",
                   control.predictor = list(compute = TRUE))
m = result.pred$marginals.linear.predictor[[101]]  
round(result.pred$summary.linear.predictor[101,], 3)

plot(m)
lines(inla.smarginal(m))