## 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
lines(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)

## 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)