knitr::opts_chunk$set(warning = FALSE, comment = "#>", fig.height = 5, fig.width = 7) library(pavpop)
In this vignette we will consider some simple examples of data that needs registration.
First, we generate a dataset with randomly shifted curves with noise
# Number of samples n <- 30 # Number of observation points m <- 30 # Observation points t <- seq(0, 1, length = m + 2)[2:(m + 1)] # Mean function theta <- function(t) dnorm(t, mean = 0.3, sd = 0.2) + dnorm(t, mean = 0.7, sd = 0.15) # Generate shifts w <- runif(n, min = -0.1, max = 0.1) # Generate data with random shifts sigma <- 0.1 y <- lapply(w, function(w) {theta(t + w) + rnorm(m, sd = sigma)}) t <- lapply(1:n, function(x) t) # Plot shifted curves plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 't', ylab = 'y(t)') legend(0.9, 3, 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]], theta(t[[1]]), lwd = 2, lty = 2)
We now set up the pavpop model to estimate in the model
# Set up basis function kts <- seq(-0.2, 1.2, length = 15)[2:14] basis_fct <- make_basis_fct(kts = kts, intercept = TRUE, control = list(boundary = c(-0.2, 1.2))) # Set up warp function warp_fct <- make_warp_fct(type = 'shift') # Estimate in the model res <- pavpop(y, t, basis_fct, warp_fct, amp_cov = NULL, warp_cov = NULL, iter = c(1, 5))
We can not plot the results using a coloring similar to before.
plot(w, res$w, xlab = 'true shifts', ylab = 'estimated shifts', pch = 19) abline(0, 1, lty = 2) t_plot <- seq(-0.2, 1.2, length = 100) plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 'warped t', ylab = expression(theta(t))) for (i in 1:n) lines(t[[i]] + res$w[i], y[[i]], col = rainbow(n)[i]) lines(t_plot, basis_fct(t_plot) %*% res$c, lwd = 2, lty = 2, col = 'red') lines(t_plot, theta(t_plot), lwd = 2, lty = 2)
We use the same data as before, but now we generate random shifts from a normal distribution
# Generate shifts scale <- 0.5 w <- rnorm(n, sd = sigma * sqrt(scale)) # Generate data with random shifts y <- lapply(1:n, function(i) {theta(t[[i]] + w[i]) + rnorm(m, sd = sigma)}) # Plot shifted curves plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 't', ylab = expression(theta(t))) for (i in 1:n) lines(t[[i]], y[[i]], col = rainbow(n)[i]) lines(t[[1]], theta(t[[1]]), lwd = 2, lty = 2)
Let us first ignore the uncertainty and use the same model as in Example 1.
res_no_uncert <- pavpop(y, t, basis_fct, warp_fct, amp_cov = NULL, warp_cov = NULL, iter = c(1, 5))
Correct specification of the model where the shifts are a normal random effect can be done as follows
warp_cov <- make_cov_fct(id_cov, noise = FALSE) res <- pavpop(y, t, basis_fct, warp_fct, amp_cov = NULL, warp_cov = warp_cov, iter = c(5, 5))
We can not plot the results, where the predicted shifts are green and the estimated shifts from the model with warps as fixed effects are red.
plot(w, res$w, xlab = 'true shifts', ylab = 'estimated/predicted shifts', pch = 19, col = 'green') points(w, res_no_uncert$w, pch = 19, col = 'red', cex = 0.5) abline(0, 1, lty = 2) t_plot <- seq(-0.2, 1.2, length = 100) plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 'warped t', ylab = expression(theta(t))) for (i in 1:n) lines(t[[i]] + res$w[i], y[[i]], col = rainbow(n)[i]) lines(t_plot, basis_fct(t_plot) %*% res$c, lwd = 2, lty = 2, col = 'red') lines(t_plot, theta(t_plot), lwd = 2, lty = 2) # Compare estimated variance parameters of warps # True standard deviation sigma * sqrt(scale) # Standard deviation of true shifts sd(w) # Estimated standard deviation res$sigma * sqrt(res$warp_cov_par) # Standard deviation of predicted warps sd(res$w)
We use the same data as in Example 2, but instad of iid normally distributed noise we now add a Matérn process to the observations.
# Generate shifts scale <- 0.5 w <- rnorm(n, sd = sigma * sqrt(scale)) amp_cov <- make_cov_fct(Matern, noise = FALSE) # Generate data with random shifts gen_dat <- function(i) { theta(t[[i]] + w[i]) + sigma * t(chol(amp_cov(t[[i]], c(100, 0.3, 2)))) %*% rnorm(m) } y <- lapply(1:n, gen_dat) # Plot shifted curves plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 't', ylab = expression(theta(t))) for (i in 1:n) lines(t[[i]], y[[i]], col = rainbow(n)[i]) lines(t[[1]], theta(t[[1]]), lwd = 2, lty = 2)
Let us first ignore the uncertainty and use the same model as in Example 1.
res_no_uncert <- pavpop(y, t, basis_fct, warp_fct, amp_cov = NULL, warp_cov = NULL, iter = c(1, 10))
And now compare to the correct specification of the model
res <- pavpop(y, t, basis_fct, warp_fct, amp_cov = amp_cov, warp_cov = warp_cov, iter = c(10, 5))
We can not plot the results
plot(w, res$w, xlab = 'true shifts', ylab = 'estimated/predicted shifts', pch = 19, col = 'green', ylim = range(res$w, res_no_uncert$w)) points(w, res_no_uncert$w, pch = 19, col = 'red', cex = 0.5) abline(0, 1, lty = 2) par(mfrow = c(1, 2)) plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 'warped t', ylab = expression(theta(t)), main = 'No random effects') for (i in 1:n) lines(t[[i]] + res_no_uncert$w[i], y[[i]], col = rainbow(n)[i]) lines(t_plot, basis_fct(t_plot) %*% res_no_uncert$c, lwd = 2, lty = 2, col = 'blue') lines(t_plot, theta(t_plot), lwd = 2, lty = 2) plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', xlab = 'warped t', ylab = expression(theta(t)), main = 'pavpop') for (i in 1:n) lines(t[[i]] + res$w[i], y[[i]], col = rainbow(n)[i]) lines(t_plot, basis_fct(t_plot) %*% res$c, lwd = 2, lty = 2, col = 'blue') lines(t_plot, theta(t_plot), lwd = 2, lty = 2) # Compare estimated variance parameters of warps # True standard deviation sigma * sqrt(scale) # Standard deviation of true shifts sd(w) # Estimated standard deviation res$sigma * sqrt(res$warp_cov_par) # Standard deviation of predicted warps sd(res$w)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.