Nothing
#' Gaussian model of functional data
#'
#' This function models the functional data using a Gaussian model extracted from
#' the principal components of the srvfs
#'
#' @param warp_data fdawarp object from [time_warping] of aligned data
#' @param n number of random samples (n = 1)
#' @param sort_samples sort samples (default = F)
#' @return Returns a fdawarp object containing \item{fs}{random aligned samples}
#' \item{gams}{random warping function samples}
#' \item{ft}{random function samples}
#' @keywords pca
#' @references Tucker, J. D., Wu, W., Srivastava, A.,
#' Generative Models for Function Data using Phase and Amplitude Separation,
#' Computational Statistics and Data Analysis (2012), 10.1016/j.csda.2012.12.001.
#' @export
#' @examples
#' out1 <- gauss_model(simu_warp, n = 10)
gauss_model <- function(warp_data,n = 1,sort_samples = FALSE){
fn <- warp_data$fn
time <- warp_data$time
qn <- warp_data$qn
gam <- warp_data$warping_functions
# Parameters
no = 3
eps = .Machine$double.eps
binsize = mean(diff(time))
M = length(time)
# compute mean and covariance in q-domain
mq_new = rowMeans(qn)
id = round(length(time)/2)
m_new = sign(fn[id,])*sqrt(abs(fn[id,])) # scaled version
mqn = c(mq_new,mean(m_new))
C = stats::cov(t(rbind(qn,m_new)))
q_s = mvtnorm::rmvnorm(n,mean=mqn,sigma=C,method="svd")
q_s = t(q_s)
end = dim(q_s)[1]
# compute the correspondence to the original function domain
fs = matrix(0,M,n)
for (k in 1:n){
fs[,k] = cumtrapzmid(time,q_s[1:(end-1),k]*abs(q_s[1:(end-1),k]),sign(q_s[end,k])*(q_s[end,k]^2),id)
}
fbar = rowMeans(fn)
fsbar = rowMeans(fs)
err = kronecker(matrix(1,1,n),fbar-fsbar)
fs = fs + err
# random warping generation
rgam = randomGamma(gam,n)
gams = matrix(0,n,M)
for (k in 1:n){
gams[k,] = invertGamma(rgam[k,])
}
gams = t(gams)
# sort functions and warpings
if (sort_samples == T){
mx = apply(fs,2, max)
out_sort = sort(mx,index.return=TRUE)
seq1 = out_sort$ix
# compute the psi-function
fy = gradient(t(rgam),binsize)
psi = fy/sqrt(abs(fy)+eps)
psi = t(psi)
ip = rep(0,n)
len = rep(0,n)
for (i in 1:n){
ip[i] = rep(1,M)%*%psi[i,]/M;
len[i] = acos(rep(1,M)%*%psi[i,]/M)
}
out_sort = sort(len,index.return=TRUE)
seq2 = out_sort$ix
# combine x-variability and y-variability
ft = matrix(0,M,n)
for (k in 1:n){
tmp = stats::approx((0:(M-1))/(M-1),fs[,seq1[k]],xout = gams[,seq2[k]])
ft[,k] = tmp$y
while (is.na(ft[,k])){
rgam2 = randomGamma(gam,1)
tmp = stats::approx((0:(M-1))/(M-1),fs[,seq1[k]],xout = invertGamma(rgam2))
ft[,k] = tmp$y
}
}
}else
{
# combine x-variability and y-variability
ft = matrix(0,M,n)
for (k in 1:n){
tmp = stats::approx((0:(M-1))/(M-1),fs[,k],xout = gams[,k])
ft[,k] = tmp$y
while (is.na(ft[,k])[1]){
rgam2 = randomGamma(gam,1)
tmp = stats::approx((0:(M-1))/(M-1),fs[,k],xout = invertGamma(rgam2))
ft[,k] = tmp$y
}
}
}
warp_data$fs = fs
warp_data$gams = rgam
warp_data$ft = ft
warp_data$qs = q_s[1:(end-1),]
warp_data$rsamps=T
return(warp_data)
}
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.