## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----sim-data, eval=F---------------------------------------------------------
# set.seed(12345)
# library(refundBayes)
# library(refund)
# library(ggplot2)
# 
# ## Dimensions
# n     <- 200
# M     <- 50
# tgrid <- seq(0, 1, length.out = M)
# 
# ## True mean function and eigenfunctions
# mu_true <- sin(pi * tgrid)
# phi1    <- sqrt(2) * sin(2 * pi * tgrid)
# phi2    <- sqrt(2) * cos(2 * pi * tgrid)
# phi_true <- rbind(phi1, phi2)            # J x M
# eigen_true     <- c(2, 0.5)
# sigma_eps_true <- 0.3
# 
# ## Simulate scores and observations
# xi_true <- cbind(rnorm(n, sd = sqrt(eigen_true[1])),
#                  rnorm(n, sd = sqrt(eigen_true[2])))
# Y_mat <- matrix(rep(mu_true, n), nrow = n, byrow = TRUE) +
#          xi_true %*% phi_true +
#          matrix(rnorm(n * M, sd = sigma_eps_true), nrow = n)
# 
# dat <- data.frame(inx = 1:n)
# dat$Y_mat <- Y_mat

## ----fit-model, eval=F--------------------------------------------------------
# fit_fpca <- refundBayes::fpca_bayes(
#   formula     = Y_mat ~ 1,
#   data        = dat,
#   spline_type = "bs",
#   spline_df   = 10,
#   niter       = 1500,
#   nwarmup     = 500,
#   nchain      = 1,
#   ncores      = 1
# )

## ----plot-model, eval=F-------------------------------------------------------
# library(ggplot2)
# p <- plot(fit_fpca)        # default prob = 0.95
# p$mu                       # posterior mean function with 95% credible band
# p$efunctions               # fixed eigenfunctions used as the FPC basis
# p$evalues                  # posterior eigenvalue SD with 95% intervals
# p$sigma                    # posterior of sigma_eps (histogram)

## ----posterior-summary, eval=F------------------------------------------------
# ## Posterior mean of the mean function on the input grid
# mu_hat <- apply(fit_fpca$mu, 2, mean)
# 
# ## Pointwise 95% credible interval for the mean function
# mu_lo  <- apply(fit_fpca$mu, 2, function(x) quantile(x, 0.025))
# mu_hi  <- apply(fit_fpca$mu, 2, function(x) quantile(x, 0.975))
# 
# ## Posterior mean of the FPC eigenvalue SDs
# lambda_hat <- apply(fit_fpca$evalues, 2, mean)
# 
# ## Posterior mean of the FPC scores (n x J matrix of point estimates)
# scores_hat <- apply(fit_fpca$scores, c(2, 3), mean)
# 
# ## Posterior mean of the residual SD
# sigma_eps_hat <- mean(fit_fpca$sigma)

## ----reconstruct, eval=F------------------------------------------------------
# Phi <- fit_fpca$efunctions       # M x J  (fixed)
# mu_hat     <- apply(fit_fpca$mu,     2,      mean)   # length M
# scores_hat <- apply(fit_fpca$scores, c(2,3), mean)   # n x J
# 
# Y_hat <- matrix(rep(mu_hat, n), nrow = n, byrow = TRUE) +
#          scores_hat %*% t(Phi)

## ----freq-comparison, eval=F--------------------------------------------------
# library(refund)
# 
# fit_freq <- refund::fpca.sc(dat$Y_mat)
# ## or: fit_freq <- refund::fpca.face(dat$Y_mat)
# 
# ## Frequentist mean function vs Bayesian posterior mean
# plot(tgrid, fit_freq$mu, type = "l", lwd = 2, col = "darkorange",
#      ylab = expression(hat(mu)(t)), xlab = "t")
# lines(tgrid, mu_hat, col = "blue", lwd = 2)
# legend("topright", legend = c("Frequentist", "Bayesian"),
#        col = c("darkorange", "blue"), lwd = 2, bty = "n")

## ----inspect-code, eval=F-----------------------------------------------------
# fpca_code_only <- refundBayes::fpca_bayes(
#   formula = Y_mat ~ 1,
#   data    = dat,
#   runStan = FALSE
# )
# cat(fpca_code_only$stancode)