Nothing
## ----set-knitr-options, cache=FALSE, echo=FALSE, warning=FALSE, message=FALSE----
library("knitr")
opts_chunk$set(message = FALSE, fig.width = 5.5)
## ----message=FALSE, warning=FALSE---------------------------------------------
library(bayesdfa)
library(ggplot2)
library(dplyr)
library(rstan)
chains = 1
iter = 10
## ----simulate-data------------------------------------------------------------
set.seed(1)
sim_dat <- sim_dfa(
num_trends = 2,
num_years = 20,
num_ts = 4
)
## -----------------------------------------------------------------------------
set.seed(1)
sim_dat$x[1,] = cumsum(rnorm(n=ncol(sim_dat$x), 0, 0.1))
sim_dat$x[2,] = cumsum(rnorm(n=ncol(sim_dat$x), 0, 1))
## ----simulate-data-plot, fig.align='center', fig.cap="Simulated data, from a model with 2 latent trends and no extremes.\\label{fig:simulate-data-plot}"----
matplot(t(sim_dat$x),
type = "l",
ylab = "Response", xlab = "Time"
)
## -----------------------------------------------------------------------------
sim_dat$pred = sim_dat$Z %*% sim_dat$x
for(i in 1:nrow(sim_dat$pred)) {
for(j in 1:ncol(sim_dat$pred)) {
sim_dat$y_sim[i,j] = sim_dat$pred[i,j] + rnorm(1,0,0.1)
}
}
## ----fit-model, message=FALSE, warning=FALSE, results='hide'------------------
f1 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, scale="zscore",
iter = iter, chains = chains, thin = 1
)
r1 <- rotate_trends(f1)
## ----fit-model-2, message=FALSE, warning=FALSE, results='hide'----------------
f2 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, scale="zscore", estimate_process_sigma = TRUE,
equal_process_sigma = FALSE,
iter = iter, chains = chains, thin = 1
)
r2 <- rotate_trends(f2)
## ----fit-model-3, message=FALSE, warning=FALSE, results='hide'----------------
f3 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, scale="center", estimate_process_sigma = TRUE,
equal_process_sigma = FALSE,
iter = iter, chains = chains, thin = 1
)
r3 <- rotate_trends(f3)
## ----fit-model-4, message=FALSE, warning=FALSE, results='hide'----------------
f4 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, scale="zscore", estimate_process_sigma = TRUE,
equal_process_sigma = TRUE,
iter = iter, chains = chains, thin = 1
)
r4 <- rotate_trends(f4)
f5 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, scale="center", estimate_process_sigma = TRUE,
equal_process_sigma = TRUE,
iter = iter, chains = chains, thin = 1
)
r5 <- rotate_trends(f5)
## -----------------------------------------------------------------------------
print(round(sim_dat$Z,2))
## -----------------------------------------------------------------------------
print(round(r1$Z_rot_mean,2))
## -----------------------------------------------------------------------------
print(round(r2$Z_rot_mean,2))
## -----------------------------------------------------------------------------
print(round(r3$Z_rot_mean,2))
## -----------------------------------------------------------------------------
print(round(r4$Z_rot_mean,2))
## -----------------------------------------------------------------------------
print(round(r5$Z_rot_mean,2))
## ----echo=FALSE---------------------------------------------------------------
m = matrix(0,5,2)
m[,1] = 1:5
m[1,2] = sum((r1$Z_rot_mean - sim_dat$Z)^2)
m[2,2] = sum((r2$Z_rot_mean - sim_dat$Z)^2)
m[3,2] = sum((r3$Z_rot_mean - sim_dat$Z)^2)
m[4,2] = sum((r4$Z_rot_mean - sim_dat$Z)^2)
m[5,2] = sum((r5$Z_rot_mean - sim_dat$Z)^2)
colnames(m) = c("Model", "RMSE-loadings")
knitr::kable(m)
## ----echo=FALSE---------------------------------------------------------------
m = matrix(0,5,2)
m[,1] = 1:5
m[1,2] = sum((r1$trends_mean - sim_dat$x)^2)
m[2,2] = sum((r2$trends_mean - sim_dat$x)^2)
m[3,2] = sum((r3$trends_mean - sim_dat$x)^2)
m[4,2] = sum((r4$trends_mean - sim_dat$x)^2)
m[5,2] = sum((r5$trends_mean - sim_dat$x)^2)
colnames(m) = c("Model", "RMSE-trends")
knitr::kable(m)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.