#' Calculates a matrix in which each column is a skewed Gaussian
#'
#' @description Calculates a matrix in which each column is a skewed Gaussian.
#' Like \code{calcEhiergaus} from TIMP package but uses a vector not a list of
#' parameter estimates.
#'
#' @param theta vector of parameter estimates
#' @param lambda wavelengths at which to calculate model
#'
#' @return \code{matrix}
#' @importFrom TIMP skew l2nu
#' @export
#'
calcE <- function(theta, lambda) {
nspec <- length(theta) / 3
nl <- length(lambda)
npare <- 3
spec <- matrix(nrow = nl, ncol = nspec)
for (i in 1:nspec) {
ioff <- (i - 1) * npare
spec[, i] <- skew(theta[ioff + 1],
theta[ioff + 2],
theta[ioff + 3],
l2nu(lambda),
nupower = 1
)
}
spec
}
#' Spectrotemporal model
#'
#' @param sim object of class \code{dat} representing data
#' @param model object of class \code{dat} representing a model
#' @param iter integer number of iterations
#' @param kroncol object of class \code{logical} that is
#' \code{TRUE} if the \code{kroneckcol} function should be used to formulate the
#' model and \code{FALSE} if the standard \code{kronecker} is to be used instead
#' @param lin defines the range to plot linearly
#' (from -\code{lin} to +\code{lin})
#' @param l_posk object of class \code{logical} indicating whether
#' positivity-constraints are enforced on the rate parameters
#'
#' @importFrom TIMP compModel
#' @importFrom stats nls nls.control
#' @importFrom utils head tail
#' @export
#'
spectemp <- function(sim,
model,
iter,
kroncol = FALSE,
lin = NA,
l_posk = FALSE) {
psisim <- as.vector(sim@psi.df)
dummy <- as.data.frame(psisim)
specpar <- unlist(model@specpar)
kinpar <- model@kinpar
seqmod <- model@seqmod
irfpar <- model@irfpar
irf <- model@irf
x2 <- sim@x2
x <- sim@x
## not using the calcEhiergaus from TIMP because it requires a
## list of specpar, as opposed to a vector
if (l_posk) {
if (kroncol) {
kronform <- psisim ~ kroneckercol(
A = calcE(sp, x2),
B = compModel(
k = exp(k),
x = x,
seqmod = seqmod,
irf = irf,
irfpar = irfpar
)
)
} else {
kronform <- psisim ~ kronecker(
calcE(sp, x2),
compModel(
k = exp(k),
x = x,
seqmod = seqmod,
irf = irf,
irfpar = irfpar
)
)
}
# NB warnOnly=TRUE, to return also in case of nonconvergence
onls <- nls(
kronform,
control = nls.control(
printEval = TRUE,
warnOnly = TRUE,
maxiter = iter
),
start = list(
k = log(kinpar),
sp = specpar,
irfpar = irfpar
),
algorithm = "plinear",
trace = TRUE
)
} else {
# not lposk
if (kroncol) {
kronform <- psisim ~ kroneckercol(
A = calcE(sp, x2),
B = compModel(
k = k,
x = x,
seqmod = seqmod,
irf = irf,
irfpar = irfpar
)
)
} else {
kronform <- psisim ~ kronecker(
calcE(sp, x2),
compModel(
k = k,
x = x,
seqmod = seqmod,
irf = irf,
irfpar = irfpar
)
)
}
# NB warnOnly=TRUE, to return also in case of nonconvergence
onls <- nls(
kronform,
control = nls.control(
printEval = TRUE,
warnOnly = TRUE,
maxiter = iter
),
start = list(
k = kinpar,
sp = specpar, irfpar = irfpar
),
algorithm = "plinear",
trace = TRUE
)
} # unexpected closing bracket
# sumonlskron <- summary(onls, correlation=TRUE) if(kroncol)
# { assign('sumonlssingle', sumonlskron, envir=.GlobalEnv) }
# else { assign('sumonlsmultiple', sumonlskron,
# envir=.GlobalEnv) }
## standardize the storage of par. est., assuming order is
## kinpar, specpar irfpar
theta <- theta()
if (length(model@irfpar) > 0) {
pars <- head(onls$m$getPars(), length(onls$m$getPars()) -
2)
theta@irfpar <- tail(onls$m$getPars(), 2)
} else {
pars <- onls$m$getPars()
}
theta@specpar <- list(tail(pars, length(specpar)))
if (l_posk) {
theta@kinpar <- exp(head(pars, length(kinpar)))
} else {
theta@kinpar <- head(pars, length(kinpar))
}
# plotterforGUI(modtype='spectemp', data=sim, model=model,
# theta=theta, result=onls, lin = lin)
list(theta = theta, onls = onls)
}
#' kroneckercol: column-wise kronecker product
#'
#' @description The column-wise kronecker product is also called
#' the Khatri–Rao product
#'
#' @param A numerical matrix
#' @param B numerical matrix
#'
#' @return column-wise kronecker product of A and B
#' @export
#'
"kroneckercol" <- function(A, B) {
ncomp <- ncol(B)
nt <- nrow(B)
nl <- nrow(A)
ar <- rep(1:nl, rep(nt, nl))
br <- rep(1:nt, nl)
ac <- 1:ncomp
A[ar, ac] * B[br, ac]
}
#' Generate linlog tics for a linear-logarithmic axis
#'
#' @param x values for which to calculate a linlog axis
#' @param mu center of axis in the original \code{x} axis
#' @param alpha linear part
#'
#' @return Returns matrix with new \code{x} values in first column and the
#' corresponding labels in the second column.
#'
"linlogtics" <- function(x, mu, alpha) {
maxorigx <- max(x)
minorigx <- min(x)
ticsl <- c(-alpha)
tics <- c(-alpha)
cntmin <- -alpha
while (cntmin > minorigx) {
cntmin <- cntmin * 10
ticsl <- append(ticsl, cntmin)
tics <- append(tics, -alpha - (alpha * log10(-cntmin / alpha)))
}
ticsl <- append(sort(ticsl), c(0, alpha))
tics <- append(sort(tics), c(0, alpha))
cntmax <- alpha
while (cntmax < maxorigx) {
cntmax <- cntmax * 10
ticsl <- append(ticsl, cntmax)
tics <- append(tics, alpha + (alpha * log10(cntmax / alpha)))
}
## new x values as column 1 new x labels as colum 2
ret <- cbind(tics, ticsl)
ret
}
#' Master plot function for paramGUI
#'
#' @param modtype either 'kin', 'spec' or 'spectemp'
#' @param X matrix of conditionally linear parameters, if any
#' @param data object of class \code{dat} containing data
#' @param model object of class \code{dat} containing data
#' @param theta object of class \code{theta} containing parameters
#' @param result object returned by \code{fitModel} or in the case
#' \code{modtype=='spectemp'}, by \code{nls}
#' @param lin The linear range for the concentration plot
#' @param mu The center of the lin-log axis is \code{lin} is specified
#' @param guessIRF Boolean to indicate whether to try and guess the
#' location of the IRF
#'
#' @return graphics
#' @importFrom fields image.plot
#' @importFrom TIMP compModel getSpecList parEst linloglines linlogplot
#' @importFrom TIMP irfparF matlinlogplot
#' @importFrom graphics abline axis barplot image lines matlines matplot
#' @importFrom graphics mtext par plot
#' @importFrom stats dnorm
#' @export
#'
"plotterforGUI" <- function(modtype = "kin",
X = matrix(),
data,
model,
theta = vector(),
result,
lin = NA,
mu = 0,
guessIRF = FALSE) {
# hard coded color table, avoiding explicit dependency on colorspace
ccs <- c(
"#A84D63",
"#AC576A",
"#B06072",
"#B46979",
"#B87280",
"#BC7A87",
"#C0838F",
"#C38B96",
"#C7939D",
"#CA9BA4",
"#CDA3AB",
"#D0ABB2",
"#D3B3B8",
"#D6BABF",
"#D8C1C5",
"#DAC8CB",
"#DDCFD1",
"#DFD5D7",
"#E0DBDC",
"#E2E0E1",
"#E0E1E0",
"#DBDDDA",
"#D4D9D3",
"#CDD4CB",
"#C5D0C3",
"#BDCBBA",
"#B5C5B1",
"#ADC0A8",
"#A4BB9E",
"#9CB594",
"#93AF8A",
"#8AA980",
"#81A475",
"#779E6A",
"#6E985F",
"#649253",
"#5A8B46",
"#508538",
"#457F27",
"#3A790D"
)
# To regenerate the color table, use the following code:
# # Note this requires:
# # @importFrom colorspace diverge_hcl
# ccs <- diverge_hcl(40,
# h = c(0, 120), c = 60, l = c(45, 90),
# power = 1.2
# )
## note that result is the return value of fitModel if
## modtype=='spec' or modtype == 'kin', and the return value
## of nls otherwise
if (!is.null(model)) {
## get clp, if modtype == 'kin' || modtype == 'spec'
if (!modtype == "spectemp") {
X <- getSpecList(result$currModel, result$currTheta)[[1]]
theta <- result$currTheta[[1]]
}
} else {
modtype <- "kin"
}
nt <- data@nt
nl <- data@nl
x <- data@x
x2 <- data@x2
observed <- data@psi.df
svdobserved <- svd(observed)
if (!is.null(model)) {
if (modtype == "kin" && length(model@irfpar) > 0) {
mu <- unlist(parEst(
result,
param = "irfpar",
dataset = 1,
verbose = FALSE
))[1]
} else {
if (modtype == "spectemp" && length(theta@irfpar) > 0) {
mu <- head(theta@irfpar, 1)
} else {
mu <- 0
}
}
} else if (guessIRF) {
lsv1 <- svd(data@psi.df)$u[, 1]
mu <- data@x[[floor((which(lsv1 == min(lsv1)) +
which(lsv1 == max(lsv1))) / 2)]]
}
op <- par(no.readonly = TRUE)
# CHANGE PLOT OPTIONS
# TODO: change plotting options to one of the alternatives:
# layout or split.screen.
if (!is.null(model)) {
if ((nt == 1) || (nl == 1)) {
par(mfrow = c(2, 2), oma = c(0, 0, 3, 0))
} else {
par(mfrow = c(3, 4), oma = c(0, 0, 3, 0))
}
} else {
par(mfrow = c(2, 3), oma = c(0, 0, 3, 0))
}
## get labels and time axis for linlog plots
dolinlog <- FALSE
if (is.na(lin)) {
xnew <- x
if (!modtype == "spec") {
if ((!is.null(model) && length(model@irfpar) > 0)) {
lin <- max(pretty(abs(mu) * 10))
} else {
lin <- max(data@x)
}
} else {
lin <- max(data@x)
}
} else {
dolinlog <- TRUE
xnew <- linloglines(x, mu = mu, alpha = lin)
newlab <- linlogtics(x, mu = 0, alpha = lin)
}
if (!is.null(model)) {
if (!modtype == "spectemp") {
residuals <- matrix(nrow = nt, ncol = nl)
residlist <- result$currModel@fit@resultlist[[1]]@resid
for (j in seq_along(residlist)) {
if (modtype == "kin") {
residuals[, j] <- residlist[[j]]
} else {
residuals[j, ] <- residlist[[j]]
}
}
} else {
residuals <- result$m$resid()
dim(residuals) <- c(nt, nl)
}
}
# PLOT DATA
if (nt == 1) {
plot(
x = x2,
y = observed,
xlab = "wavelength (nm)",
ylab = "",
main = "Data",
type = "l",
xlim = c(
min(x2),
max(x2)
)
)
if (!is.null(model)) {
lines(
x = x2,
y = observed - residuals[1, ],
col = "red"
)
}
abline(0, 0)
} else if (nl == 1) {
## linlogplot(x=x, y=observed, mu=mu, alpha=lin, xlab = 'time
## (ps)',ylab='', main = 'Data', type = 'l', xlim=c(min(x),
## max(x)))
plot(
x = x,
y = observed,
xlab = "time (ps)",
ylab = "",
main = "Data",
type = "l",
xlim = c(min(x), max(x))
)
if (!is.null(model)) {
lines(
x = x,
y = observed - residuals[, 1],
col = "red"
)
}
abline(0, 0)
} else {
m <- par("mar")
par(mar = c(m[1:3], 3))
if (dolinlog) {
image(
xnew,
x2,
observed,
ylab = "wavelength (nm)",
xaxt = "n",
main = "Data",
xlab = "time (ps)",
col = ccs
)
axis(1, at = newlab[, 1], labels = newlab[, 2])
# mtext(side = 1, newlab[,2], at= newlab[,1], line = 1)
image.plot(xnew,
x2,
observed,
legend.only = TRUE,
col = ccs
)
} else {
image.plot(
xnew,
x2,
observed,
ylab = "wavelength (nm)",
main = "Data",
xlab = "time (ps)",
col = ccs
)
}
par(mar = m)
# PLOT SVD DATA
plot(log10(svdobserved$d),
main = "log(sing val. data)",
ylab = ""
)
lsv1 <- svdobserved$u[, 1]
# maxlsv1<-max(lsv1) minlsv1<-min(lsv1)
# extrlsv1<-max(abs(maxlsv1),abs(minlsv1)) if (minlsv1 > 0)
# minlsv1=0 if (maxlsv1 < 0) maxlsv1=0
if (dolinlog) {
linlogplot(
x = x,
y = lsv1,
mu = mu,
alpha = lin,
main = "1st LSV data",
xlab = "time (ps)",
type = "l",
ylab = "",
xlim = c(min(x), max(x)),
ylim = c(min(
lsv1,
0
), max(lsv1, 0))
)
} else {
plot(
x = x,
y = lsv1,
main = "1st LSV data",
xlab = "time (ps)",
type = "l",
ylab = "",
xlim = c(min(x), max(x)),
ylim = c(min(lsv1, 0), max(lsv1, 0))
)
}
abline(0, 0, lty = 3)
plot(
x2,
svdobserved$v[, 1],
main = "1st RSV data",
xlab = "wavelength (nm)",
type = "l",
ylab = ""
)
abline(0, 0, lty = 3)
# moved plotting of 2nd LSV and RSV after plotting of spectra
}
if (!is.null(model)) {
if (nt > 1) {
# PLOT CONCENTRATIONS
if (!modtype == "spec") {
C <- compModel(
k = theta@kinpar,
x = x,
irfpar = theta@irfpar,
irf = model@irf,
seqmod = model@seqmod
)
} else {
# modtype == 'spec'
C <- X
# assign('Cest',cbind(x, C), envir = globalenv())
}
# if(lin == 0) # reversed plots
if (length(theta@irfpar) > 0) {
if (theta@irfpar[2] > 0) {
irf_for_plotting <- dnorm(x, theta@irfpar[1], theta@irfpar[2])
irf_for_plotting <-
irf_for_plotting / max(irf_for_plotting) * max(C)
} else {
irf_for_plotting <- rep(0, length(x))
}
}
if (dolinlog) {
matlinlogplot(
x,
C,
mu,
lin,
ylab = "",
xlab = "time (ps)",
main = "Concentrations",
type = "l",
lty = 1
)
if (length(theta@irfpar) > 0) {
matlinlogplot(
x,
irf_for_plotting,
mu,
lin,
type = "l",
lty = 2,
add = TRUE
)
}
} else {
matplot(
x,
C,
xlab = "time (ps)",
ylab = "",
main = "Concentrations",
type = "l",
lty = 1
)
if (length(theta@irfpar) > 0) {
matplot(
x,
irf_for_plotting,
type = "l",
lty = 2,
add = TRUE
)
}
}
if (data@simdata) {
if ((modtype == "kin") || (modtype == "spectemp")) {
amplitudes <- data@amplitudes
ncomp <- length(amplitudes)
aC2 <- data@C2 %*% diag(
1 / amplitudes, ncomp,
ncomp
)
if (dolinlog) {
matlinlogplot(
x = x,
mu = mu,
alpha = lin,
y = aC2,
col = "blue",
lty = 3,
add = TRUE,
type = "l"
)
} else {
matplot(
x = x,
y = aC2,
col = "blue",
lty = 3,
add = TRUE,
type = "l"
)
}
} else {
if (dolinlog) {
matlinlogplot(
x = x,
mu = mu,
alpha = lin,
y = data@C2,
col = "blue",
lty = 3,
add = TRUE,
type = "l"
)
} else {
matplot(
x = x,
y = data@C2,
col = "blue",
lty = 3,
add = TRUE,
type = "l"
)
}
}
}
abline(0, 0, lty = 3)
} else {
barplot(
X[1, ],
main = "Amplitudes",
ylab = "",
xlab = "component",
lty = 1
)
}
if (nl > 1) {
# PLOT SPECTRA
if (modtype == "kin") {
E <- X
# assign('Eest',cbind(x2, E), envir = globalenv())
}
if (modtype == "spec") {
E <- calcEhiergaus(
lambda = x2,
theta = theta@specpar,
nupower = 1
)
}
if (modtype == "spectemp") {
E <- calcE(lambda = x2, theta = unlist(theta@specpar))
}
matplot(
x2,
E,
main = "Spectra",
ylab = "",
xlab = "wavelength (nm)",
type = "l",
lty = 1
)
abline(0, 0, lty = 3)
if (data@simdata) {
if (modtype == "kin") {
aE2 <- data@E2 %*% diag(
amplitudes, ncomp,
ncomp
)
matlines(x2, aE2, lty = 3, col = "blue")
} else {
matlines(x2, data@E2, lty = 3, col = "blue")
}
}
} else {
barplot(
X[1, ],
main = "Amplitudes",
ylab = "",
xlab = "component",
lty = 1
)
}
# Moved 2nd LSV and RSV here:
if (nt > 1 && nl > 1) {
if (dolinlog) {
linlogplot(
x = x,
y = svdobserved$u[, 2],
mu = mu,
alpha = lin,
main = "2nd LSV data",
xlab = "time (ps)",
type = "l",
ylab = "",
xlim = c(min(x), max(x))
)
} else {
plot(
x = x,
y = svdobserved$u[, 2],
main = "2nd LSV data",
xlab = "time (ps)",
type = "l",
ylab = "",
xlim = c(
min(x),
max(x)
)
)
}
abline(0, 0, lty = 3)
plot(
x2,
svdobserved$v[, 2],
main = "2nd RSV data",
xlab = "wavelength (nm)",
type = "l",
ylab = ""
)
abline(0, 0, lty = 3)
}
# PLOT RESIDS
if (nt == 1) {
plot(
x = x2,
y = residuals[1, ],
xlab = "wavelength (nm)",
ylab = "",
main = "Residuals",
type = "l",
xlim = c(
min(x2),
max(x2)
)
)
abline(0, 0)
} else {
if (nl == 1) {
linlogplot(
x = x,
y = residuals[, 1],
mu = mu,
alpha = lin,
xlab = "time (ps)",
ylab = "",
main = "Residuals",
type = "l",
xlim = c(
min(x),
max(x)
)
)
abline(0, 0)
} else {
m <- par("mar")
par(mar = c(m[1:3], 3))
if (dolinlog) {
image(
xnew,
x2,
residuals,
ylab = "wavelength (nm)",
xaxt = "n",
main = "Residuals",
col = ccs,
xlab = "time (ps)"
)
axis(1, at = newlab[, 1], labels = newlab[
,
2
])
image.plot(
xnew,
x2,
residuals,
ylab = "wavelength (nm)",
legend.only = TRUE,
main = "Residuals",
col = ccs
)
} else {
image.plot(
xnew,
x2,
residuals,
ylab = "wavelength (nm)",
main = "Residuals",
col = ccs,
xlab = "time (ps)"
)
}
par(mar = m)
svdresid <- svd(residuals)
logsvd <- log10(svdresid$d)
lenlogsvd <- length(logsvd)
nsingvalres <- lenlogsvd - dim(E)[2]
plot(logsvd[1:nsingvalres],
main = "log(sing. val. resid.)",
ylab = ""
)
linlogplot(
x = x,
y = svdresid$u[, 1],
mu = mu,
alpha = lin,
main = "1st LSV resid.",
xlab = "time (ps)",
type = "l",
ylab = "",
xlim = c(min(x), max(x))
)
abline(0, 0, lty = 3)
plot(
x2,
svdresid$v[, 1],
main = "1st RSV resid",
xlab = "wavelength (nm)",
type = "l",
ylab = ""
)
abline(0, 0, lty = 3)
}
}
# ADD TITLE
if (!modtype == "spec") {
kinest <- paste("Kin par:", toString(signif(theta@kinpar,
digits = 4
)))
if (length(theta@irfpar) > 0) {
kinest <- paste(kinest, " IRF par:", toString(signif(theta@irfpar,
digits = 4
)))
}
mtext(kinest,
side = 3,
outer = TRUE,
line = 1
)
}
if (!modtype == "kin") {
specest <- paste("Spec par:", toString(signif(
unlist(
theta@specpar
),
digits = 4
)))
mtext(specest,
side = 3,
outer = TRUE,
line = -0.5
)
}
}
# RESET FORMER PLOT OPTIONS
par(op)
# try 20111104 tkdestroy(winmodel) this only kills one
# winmodel ...
}
#' Simulate data
#'
#' @description Calculates an object of class 'kin'. <TODO>
#'
#' @author Katharine M. Mullen
#' @author Ivo H. M. van Stokkum
#'
#' @param kinpar vector of rate constants
#' @param tmax last time point
#' @param deltat time step
#' @param specpar vector of spectral parameters for location, width, skewness
#' @param lmin minimum wavelength (nm)
#' @param lmax maximum wavelength (nm)
#' @param deltal wavelength step
#' @param sigma noise level
#' @param irf logical for IRF usage
#' @param irfpar vector of IRF parameters for location, width
#' @param seqmod logical for sequential model
#' @param dispmu logical for dispersion of IRF location mu
#' @param nocolsums logical for <TODO>
#' @param disptau logical for dispersion of IRF width tau
#' @param parmu vector of dispersion parameters for IRF location mu
#' @param partau vector of dispersion parameters for IRF width tau
#' @param lambdac center wavelength for dispersion
#' @param fullk logical for full K matrix
#' @param kmat K matrix
#' @param jvec input vector
#' @param specfun function for spectral shape
#' @param nupow power of nu in spectral model
#' @param irffun function for IRF
#' @param kinscal vector of kinetic scaling parameters
#' @param lightregimespec <TODO>
#' @param specdisp logical for dispersion parameters of spectral parameters
#' @param specdisppar vector of dispersion parameters of spectral parameters
#' @param parmufunc <TODO>
#' @param specdispindex <TODO>
#' @param amplitudes amplitudes of components
#' @param specref <TODO>
#' @param nosiminfo logical for hiding simulation information
#'
#' @return an object of class 'kin'
#' @importFrom TIMP dat kin compModel specparF calcEhiergaus
#' @importFrom stats rnorm
#' @export
#'
"simndecay_gen_paramGUI" <-
function(kinpar,
tmax,
deltat,
specpar = vector(),
lmin,
lmax,
deltal,
sigma,
irf = FALSE,
irfpar = vector(),
seqmod = FALSE,
dispmu = FALSE,
nocolsums = FALSE,
disptau = FALSE,
parmu = list(),
partau = vector(),
lambdac = 0,
fullk = FALSE,
kmat = matrix(),
jvec = vector(),
specfun = "gaus",
nupow = 1,
irffun = "gaus",
kinscal = vector(),
lightregimespec = list(),
specdisp = FALSE,
specdisppar = list(),
parmufunc = "exp",
specdispindex = list(),
amplitudes = vector(),
specref = 0,
nosiminfo = TRUE) {
if (tmax > 0) {
x <- seq(0, tmax, deltat)
} else {
x <- c(1)
}
nt <- length(x)
ncomp <- length(kinpar)
x2 <- seq(lmin, lmax, deltal)
nl <- length(x2)
if (specdisp) {
## store all the spectra; could do it otherwise if mem. is an
## issue
EList <- list()
for (i in 1:nt) {
sp <- specparF(
specpar = specpar,
xi = x[i],
i = i,
specref = specref,
specdispindex = specdispindex,
specdisppar = specdisppar,
parmufunc = parmufunc
)
EList[[i]] <- calcEhiergaus(sp, x2, nupow)
}
} else if (lmin == lmax) {
E2 <- matrix(1, nrow = 1, ncol = ncomp)
# TODO: set modType to 0?
} else {
E2 <- calcEhiergaus(specpar, x2, nupow)
}
if (!(dispmu || disptau)) {
if (nt == 1) {
C2 <- matrix(amplitudes, nrow = 1, ncol = ncomp)
# TODO: set modType to 0?
} else {
C2 <- compModel(
k = kinpar,
x = x,
irfpar = irfpar,
irf = irf,
seqmod = seqmod,
fullk = fullk,
kmat = kmat,
jvec = jvec,
amplitudes = amplitudes,
lightregimespec = lightregimespec,
nocolsums = nocolsums,
kinscal = kinscal
)
}
if (specdisp) {
psisim <- matrix(nrow = nt, ncol = nl)
E2 <- EList[[1]]
for (i in 1:nt) {
psisim[i, ] <- t(as.matrix(C2[i, ])) %*% t(EList[[i]])
}
} else {
psisim <- C2 %*% t(E2)
}
} else {
psisim <- matrix(nrow = nt, ncol = nl)
for (i in 1:nl) {
irfvec <- irfparF(
irfpar,
lambdac,
x2[i],
i,
dispmu,
parmu,
disptau,
partau,
"",
"",
"gaus"
)
C2 <- compModel(
k = kinpar,
x = x,
irfpar = irfpar,
irf = irf,
seqmod = seqmod,
fullk = fullk,
kmat = kmat,
jvec = jvec,
amplitudes = amplitudes,
lightregimespec = lightregimespec,
nocolsums = nocolsums,
kinscal = kinscal
)
psisim[, i] <- C2 %*% cbind(E2[i, ])
}
}
dim(psisim) <- c(nt * nl, 1)
psi.df <- psisim + sigma * rnorm(nt * nl)
dim(psi.df) <- c(nt, nl)
if (nosiminfo) {
dat(
psi.df = psi.df,
x = x,
nt = nt,
x2 = x2,
nl = nl,
simdata = FALSE
)
} else {
kin(
psi.df = psi.df,
x = x,
nt = nt,
x2 = x2,
nl = nl,
C2 = C2,
E2 = E2,
kinpar = kinpar,
specpar = specpar,
seqmod = seqmod,
irf = irf,
irfpar = irfpar,
dispmu = dispmu,
disptau = disptau,
parmu = parmu,
partau = partau,
lambdac = lambdac,
simdata = TRUE,
fullk = fullk,
kmat = kmat,
jvec = jvec,
amplitudes = amplitudes
)
}
}
# '.onLoad' <- function (lib, pkg) { startGUI() }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.