knitr::opts_chunk$set(warning = FALSE, comment = "#>", fig.height = 5, fig.width = 7) library(pavpop)
# Number of samples n <- 20 # Number of observation points m <- 100 # Observation points t <- seq(0, 1, length = m + 2)[2:(m + 1)] # Common basis function (both mean and amplitude variation) kts <- seq(0, 1, length = 10) basis_fct <- make_basis_fct(kts = kts, intercept = TRUE, control = list(boundary = c(-0.5, 1.5))) amp_fct <- make_basis_fct(type = 'intercept') df <- attr(basis_fct, 'df') # Generate true mean weights beta_t <- rexp(df, 0.5) # Generate random intercepts b_t <- rnorm(n, sd = 2) # Generate warping function and random parameters sigma <- 0.1 warp_fct <- make_warp_fct(type = 'linear') warp_cov_t <- matrix(c(2, 0.5, 0.5, 1), 2, 2) w_t <- replicate(n, (t(chol(warp_cov_t)) %*% rnorm(2, sd = sigma))[, 1]) # Generate data y <- lapply(1:n, function(i) {as.numeric(basis_fct(warp_fct(w_t[, i], t)) %*% beta_t + amp_fct(t) %*% b_t[i] + rnorm(m, sd = sigma))}) t <- lapply(1:n, function(x) t) # Plot observations plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 't', ylab = 'y(t)') legend(0.7, range(y)[2], legend = expression(theta(t)), lty = 2, lwd = 2) for (i in 1:n) lines(t[[i]], y[[i]], col = rainbow(n)[i]) lines(t[[1]], basis_fct(t[[1]]) %*% beta_t , lwd = 2, lty = 2)
We now set up pavpop to estimate in the sitar model
amp_cov <- make_cov_fct(id_cov, noise = FALSE) warp_cov <- make_cov_fct(unstr_cov, param = c(1, 1, 0), noise = FALSE) res <- pavpop(y, t, basis_fct, warp_fct, amp_cov, warp_cov, amp_fct, iter = c(10, 10)) # Plot aligned samples t_p <- seq(-0.2, 1.2, length = m + 2) plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 'warped t', ylab = 'y(t)') legend(1, range(y)[2], legend = c(expression(theta(t)), expression(hat(theta)(t))), lty = 2, lwd = 2, col = c('black', 'red')) for (i in 1:n) lines(warp_fct(res$w[, i], t[[i]]), y[[i]], col = rainbow(n)[i]) lines(t_p, basis_fct(t_p) %*% beta_t , lwd = 2, lty = 2) lines(t_p, basis_fct(t_p) %*% res$c, lwd = 2, lty = 2, col = 'red') plot(as.numeric(w_t), as.numeric(res$w), xlab = 'true warp variables', ylab = 'predicted warp variables', pch = 19, col = c('green', 'blue'), asp = 1) abline(0, 1, lty = 2)
Are the parameter estimates okay?
# Noise standard deviation: 0.1 res$sigma # Amplitude standard deviations: 2 res$sigma * sqrt(res$amp_cov_par) # Warp covariance # True sigma^2 * warp_cov_t # Observed cov(t(w_t)) # estimated res$sigma^2 * warp_cov(1:2, res$warp_cov_par)
library(sitar) dat <- data.frame(x = unlist(t), y = unlist(y), id = rep(1:n, each = m)) res_sitar <- sitar(x = x, y = y, id = id, data = dat, df = 14) plot(res_sitar, opt = 'd', xlim = c(-0.2, 1.2), ylim = range(y), lwd = 2, lty = 2, col = 'red', xlab = 'warped t', ylab = 'y(t)') lines(t_p, basis_fct(t_p) %*% beta_t , lwd = 2, lty = 2) legend(1, range(y)[2], legend = c(expression(theta(t)), expression(hat(theta)(t))), lty = 2, lwd = 2, col = c('black', 'red')) # Plot warped curves abc <- random.effects(res_sitar) abc$a <- 0 warped <- with(dat, xyadj(x, y, id, res_sitar, abc)) for (i in 1:n) with(warped, lines(x[1:m + (i - 1) * m], y[1:m + (i - 1) * m], col = rainbow(n)[i])) plot(as.numeric(w_t), rbind(-abc$b * exp(abc$c), exp(abc$c) - 1), xlab = 'true warp variables', ylab = 'predicted warp variables', pch = 19, col = c('green', 'blue'), asp = 1) abline(0, 1, lty = 2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.