## Do not edit this file manually.
## It has been automatically generated from mixAR.org.
# setGeneric("coef")
# setMethod("coef", "mle", function(object) object@fullcoef )
# setMethod("coef", "summary.mle", function(object) object@coef )
setGeneric("lik_params",
function(model)
standardGeneric("lik_params"),
useAsDefault=FALSE)
setGeneric("make_fcond_lik",
function(model, ts)
standardGeneric("make_fcond_lik"),
useAsDefault=FALSE)
setGeneric("mix_hatk",
function(model, x, index, xcond)
standardGeneric("mix_hatk"),
useAsDefault=FALSE)
setGeneric("mix_ek",
function(model, x, index, xcond, scale)
standardGeneric("mix_ek"),
useAsDefault=FALSE)
setGeneric("mix_pdf",
function(model, x, index, xcond, scale)
standardGeneric("mix_pdf"),
useAsDefault=FALSE)
setGeneric("mix_cdf",
function(model, x, index, xcond, scale)
standardGeneric("mix_cdf"),
useAsDefault=FALSE)
setGeneric("mix_qf", # 2021-04-30 new
function(model, p, x, index, xcond)
standardGeneric("mix_qf"),
useAsDefault=FALSE)
setClass("MixAR",
representation(prob = "numeric",
order = "numeric",
shift = "numeric",
scale = "numeric",
arcoef = "raggedCoef"
# todo: dist, pdf, cdf, qdf
# (This note is from 2012, is it still relevant?)
),
contains="VIRTUAL"
)
## zasega default stoynosti na prob sa NA, na shift 0, na scale 1.
## todo: dali vsichki da sa NA?
setMethod("initialize",
signature(.Object = "MixAR"),
function (.Object, arcoef, order, prob, shift, scale, model, ...)#2012-10-30 new arg model
{
mis_order <- missing(order)
if(missing(model)){ # 2012-10-30 lazy; but do not want to risk breaking this now.
if(missing(arcoef)) # new 2011-12-02
arcoef <- raggedCoef(p = order) # one of 'arcoef' and 'order' must be present
if(is.list(arcoef) && any(names(arcoef) == "s")){
## 13-09-2018 Updated to create "raggedCoefS" objects
## Note: previous code not touched
if(length(arcoef) != 3)
stop("Incorrect number of elements in the list for mixSAR model.")
if(all(names(arcoef) != "a") && all(names(arcoef) != "as")){
names(arcoef)[which(names(arcoef) != "s")] <- c("a", "as")
}else{
if(all(names(arcoef) != "a"))
names(arcoef)[ which(names(arcoef) == "")] <- "a"
if(all(names(arcoef) != "as"))
names(arcoef)[ which(names(arcoef) == "")] <- "as"
}
if(is.matrix(arcoef$a)){
arcoef$a <- split(arcoef$a, row(arcoef$a))
}
if(is.matrix(arcoef$as)){
arcoef$as <- split(arcoef$as, row(arcoef$as))
}
if(length(arcoef$a) != length(arcoef$as))
stop("Incorrect imputation of coefficient")
arcoef <- new("raggedCoefS", a = arcoef$a,
as = arcoef$as, s = arcoef$s)
}else{
if(!is(arcoef, "raggedCoef")){ # 2012-09-04 dobavyam "matrix"
if(is.matrix(arcoef)){
arcoef <- split(arcoef, row(arcoef)) # turn into a list of rows
if(!mis_order){
arcoef <- lapply(1:length(arcoef),
function(x) arcoef[[x]][ seq_len(order[x]) ])
}
}
arcoef <- new("raggedCoef", arcoef) # hoping 'new' knows how to do this...
}
}
wrkorder <- row_lengths(arcoef)
if(!missing(order) &&
(length(order) != length(wrkorder) || !all(order == wrkorder))){
message("Arg's 'arcoef' and 'order' are not consistent, ignoring 'order'.")
}
order <- wrkorder
## adjust 'prob' if necessary
if(missing(prob))
prob <- rep(as.numeric(NA), length(order))
else if(length(prob) == 1)
## 2020-03-25 TODO: shouldn't this be set to equal prob's?
prob <- rep(prob, length(order))
else if(length(prob) != length(order))
stop("length of 'prob' should be 1 or the length of 'order'")
## adjust 'shift' if necessary
if(missing(shift))
shift <- rep(0, length(order))
else if(length(shift) == 1)
shift <- rep(shift, length(order))
else if(length(shift) != length(order))
stop("length of 'shift' should be 1 or the length of 'order'")
## adjust 'scale' if necessary
if(missing(scale))
scale <- rep(1, length(order))
else if(length(scale) == 1)
scale <- rep(scale, length(order))
else if(length(scale) != length(order))
stop("length of 'scale' should be 1 or the length of 'order'")
}else{ # arg. 'model' supplied
ncomp <- length(model@order)
if(missing(order))
order <- model@order
else if(length(order) != ncomp || any(order != model@order))
stop("Arg. 'order' (if present) must be the same as 'model@order'.")
## else order is same as model@order, nothing to do
if(missing(prob))
prob <- model@prob
else if(length(prob) == 1)
prob <- rep(prob, ncomp)
else if(length(prob) != ncomp)
stop("Arg. 'prob' (if present) must have the same length as 'model@prob'.")
# else use 'prob' as is
if(missing(shift))
shift <- model@shift
else if(length(shift) == 1)
shift <- rep(shift, ncomp)
else if(length(shift) != ncomp)
stop("Arg. 'shift' (if present) must have the same length as 'model@shift'.")
# else use 'shift' as is
if(missing(scale))
scale <- model@scale
else if(length(scale) == 1)
scale <- rep(scale, ncomp)
else if(length(scale) != ncomp)
stop("Arg. 'scale' (if present) must have the same length as 'model@scale'.")
# else use 'scale' as is
if(missing(arcoef))
arcoef <- model@arcoef
## else{ ##arcoef supplied on top of model
##13-09-2018 Updated to create "raggedCoefS" objects
##Note: previous code not touched
## if(is.list(arcoef) & any(names(arcoef) == "s")){
## if(length(arcoef)!=3)stop(
## "Incorrect number of elements in the list to build mixSAR model.
## See help page for details")
## if(all(names(arcoef)!="a") && all(names(arcoef)!="as")){
## names(arcoef)[which(names(arcoef)!="s")] <- c("a", "as")
## }else{
## if(all(names(arcoef)!="a")) names(arcoef)[
## which(names(arcoef) == "")] <- "a"
## if(all(names(arcoef)!="as")) names(arcoef)[
## which(names(arcoef) == "")] <- "as"
## }
## if(is.matrix(arcoef$a)){
## arcoef$a <- split(arcoef$a, row(arcoef$a))
## }
## if(is.matrix(arcoef$as)){
## arcoef$as <- split(arcoef$as, row(arcoef$as))
## }
## if(length(ar$a)!=length(ar$as))stop(
## "Incorrect imputation of coefficient")
## TODO: This is wrong, I don't see variables a and s
## arcoef <- new("raggedCoefS", a=a, as=as, s=s)
##If re-established, remember to add } before else
else{
if(!is(arcoef, "raggedCoef")){ # 2012-09-04 dobavyam "matrix"
if(is.matrix(arcoef)){
arcoef <- split(arcoef, row(arcoef)) # turn into a list of rows
if(!mis_order){
arcoef <- lapply(1:length(arcoef),
function(x) arcoef[[x]][ seq_len(order[x]) ])
}
}
arcoef <- new("raggedCoef", arcoef) # hoping 'new' knows how to do this
}
}
wrkorder <- row_lengths(arcoef)
if(length(order) != length(wrkorder) || any(order != wrkorder))
stop("Arg. 'arcoef' (if present) must be consistent with 'model@order'.")
}
callNextMethod(.Object, arcoef = arcoef, order = order,
prob = prob, shift = shift, scale = scale, ...)
}
)
# 2012-10-31 commenting out. It turns out that setClass now returns a generator which is
# naturally to be named after the class. Better to avoid confusion.
# (at the same time, assigning the return value of setClass to check)
# 2012-10-31 removing the comments. mixAR is virtual, so cannot create objects from it.
# But it makes sense to use the name for something more useful.
# this has diferent syntax from the "new"
# function. Also, the defaults are NA, not zero.
mixAR <- function(template, coef, ..., filler = NA_real_){ # todo: make more versatile
if(missing(coef))
coef <- list()
if(!missing(template)){
coef$order <- template
}
coef <- .canonic_coef(coef, filler)
g <- length(coef$order)
if(g == 1)
stop("A mixAR model must have two or more components.")
new("MixARGaussian", arcoef = coef$arcoef, order = coef$order,
prob = coef$prob, shift = coef$shift, scale = coef$scale, ...)
}
setGeneric("mixAR") # note: the name of the 1st arg. is appropriate for the default only.
setMethod("mixAR", signature(template = "MixAR"),
function(template, coef, ..., filler = NA){
model <- template
if(missing(coef))
return(model)
## TODO: need to check that number of components in coef are consistent with
## ncomp
ncomp <- .nmix(model)
for(s in c("prob", "shift", "scale", "arcoef", "order")){
val <- coef[[s]]
## TODO: more care needed here
if(!is.null(val)){
if(length(val) == 1)
val <- rep(val, ncomp)
slot(model, s) <- val
}
}
model
})
.canonic_coef <- function(coef, filler){
wrk <- coef
len <- sapply(c("prob", "shift", "scale", "arcoef", "order"),
function(x){ if(is.null(wrk[[x]]))
wrk[[x]] <<- filler # non-local assignment!
length(wrk[[x]])})
len <- unique(len)
g <- len[len > 1]
if(length(g) > 1)
stop("Conflicting lengths of arguments.")
## 2018-11-17: bug-fix: changed "<-" to "<<-"
lapply(c("prob", "shift", "scale"),
function(x) if(length(wrk[[x]]) == 1) wrk[[x]] <<- rep(wrk[[x]], g))
wrk$arcoef <- if(is.null(wrk$order))
raggedCoef(value = wrk$arcoef)
else
raggedCoef(p = wrk$order, value = wrk$arcoef)
wrk
}
setMethod("show", "MixAR", # 12-09-2018 "show" method adapted to handle class "raggedCoefS"
function(object) { #note: previous content has not been changed, only "if" were added
cl <- class(object)
g <- length(object@prob)
p <- max(object@order)
mcoefs <- NULL
if(inherits(object@arcoef, "raggedCoefS")){
ps <- max(object@arcoef@ps)
s <- object@arcoef@s
}
## cat("(To see the internal structure of the object, use function 'str'.)\n\n")
cat("An object of class \"", cl, "\"\n", sep="")
cat("Number of components:", length(object@prob), "\n")
if(p > 0){
if(inherits(object@arcoef, "raggedCoefS")){
mcoefs <- ragged2char(object@arcoef@as)
colnames(mcoefs) <- paste("ar_", seq_len(ps)*s, sep="")
wrk <- lapply(1:g, function(i) c(object@prob[i], object@shift[i],
object@scale[i], object@order[i],
object@arcoef@a[[i]]))#), object@arcoef@ps[i]))
mcoef <- cbind(ragged2char(wrk), as.character(object@arcoef@ps))
colnames(mcoef) <- c("prob", "shift", "scale", "order",
paste("ar_", seq_len(p), sep=""), "s_order")
}else{
wrk <- lapply(1:g, function(i) c(object@prob[i], object@shift[i],
object@scale[i], object@order[i],
object@arcoef@a[[i]]))
mcoef <- ragged2char(wrk)
colnames(mcoef) <- c("prob", "shift", "scale", "order",
paste("ar_", seq_len(p), sep=""))
}
rownames(mcoef) <- paste("Comp_", 1:g, sep="")
print(cbind(mcoef, mcoefs), na.print="", quote=FALSE, justify="right")
}else{
m <- cbind( prob = object@prob
, shift = object@shift
, scale = object@scale
, order = object@order
)
rownames(m) <- paste("Comp_", 1:g, sep="")
cat("The AR orders of all components are 0.\n")
print(m)
}
invisible(object)
})
setAs("MixAR", "list",
function(from){
# nams <- slotNames(from)
# structure(lapply(nams, function(x) slot(from, x)),
# names = nams)
list(prob = from@prob,
order = from@order,
shift = from@shift,
scale = from@scale,
arcoef = from@arcoef@a ) # notice @a
})
setMethod("lik_params", c(model = "MixAR"),
function(model){
k <- .nmix(model) # note: this is for estimation purposes,
# only k-1 of the mixing probs are passed on.
c( model@prob[-k], model@scale, model@shift, ragged2vec(model@arcoef)
)
})
setMethod("mix_ncomp", signature(x="MixAR"),
function(x){
length(x@prob)
})
setMethod("row_lengths", signature(x="MixAR"),
function(x){
x@order
})
## 2012-11-11 new
.nmix <- function(model) length(model@prob) # number of components in a mixAR model
setMethod("make_fcond_lik", signature(model="MixAR", ts="numeric"),
function(model, ts){ ## mixar order is assumed constant in this function.
locts <- ts
mixar <- model
k <- mix_ncomp(mixar) # 2011-07-13: was length(mixar@prob)
## n <- length(lik_params(mixar))
## todo: is it worth making this function sufficiently general, so that
## subclasses of mixAR do not need to provide their own version when they
## have additional parameters?
ind_prob <- seq_len(k-1) # these need to be in sync with lik_params
ind_scale <- k:(k+k-1)
ind_shift <- (2*k):(2*k+k-1)
ind_arcoef <- - c( ind_prob, ind_scale, ind_shift) ## note: negative index
f <- function(x){
mixar@prob <- c(x[ind_prob], 1 - sum(x[ind_prob]))
mixar@scale <- x[ind_scale]
mixar@shift <- x[ind_shift]
mixar@arcoef <- rag_modify(mixar@arcoef, x[ind_arcoef])
res <- cond_loglik(mixar, locts)
## 2020-06-12 was: cat("x = ", x, "\n")
## cat("f(x) = ", res, "\n\n")
res
}
f
})
setMethod("mix_hatk", signature(model="MixAR", x="numeric", index="numeric", xcond="missing"),
function(model, x, index, xcond) {
mixFilter(x, model@arcoef, index, shift=model@shift)
})
setMethod("mix_ek", signature(model="MixAR", x="numeric", index="numeric", xcond="missing",
scale="missing"),
function(model, x, index, xcond, scale) {
mixFilter(x, model@arcoef, index, shift=model@shift,
residual=TRUE)
})
setMethod("mix_ek", signature(model="MixAR", x="numeric", index="numeric", xcond="missing",
scale="logical"),
function(model, x, index, xcond, scale) {
mixFilter(x, model@arcoef, index, shift=model@shift,
residual=TRUE, scale=model@scale)
})
## kogato xcond e dadeno, x se interpretira kato stoynosti na x_t za koito da se smyata
## Prognoz se smyata samo za xcond. Zasega samo edna stoynost na xcond.
## Primerno prilozhenie - plot na pdf ili cdf pri dadeni xcond.
setMethod("mix_ek", signature(model="MixAR", x="numeric", index="missing", xcond="numeric",
scale="missing"),
function(model, x, index, xcond) {
wrk <- mix_hatk(model, xcond, length(xcond)+1)
x - wrk ## todo: check!
})
setMethod("mix_ek", signature(model="MixAR", x="numeric", index="missing", xcond="numeric",
scale="logical"),
function(model, x, index, xcond, scale) {
wrk <- mix_hatk(model, xcond, length(xcond)+1)
(x - wrk) / model@scale
})
mixARgen <- setClass("MixARgen", representation(dist = "list" ), contains="MixAR")
setMethod("initialize",
signature(.Object = "MixARgen"),
function (.Object, dist, model, ...) # 2012-10-30 new arg. 'model'
{ # todo: investigate the following:
# .Object <- callNextMethod() # predava samo '...' args, seemingly contradicts doc.
.Object <- if(missing(model)) callNextMethod()
else callNextMethod(.Object, model = model, ...)
if(!missing(dist)){
if(is.list(dist))
.Object@dist <- dist
else
stop("Argument 'dist' must be a list.")
}else if(.hasSlot(model, "dist"))
.Object@dist <- model@dist
# else if(is(model, "MixARGaussian")) TODO: makes sense to set Gaussian components.
else
stop("Cannot set the distribution since argument 'dist' is missing.")
.Object
}
)
setMethod("show", "MixARgen", # todo: this is only a quick fix, more work is needed
function(object) {
callNextMethod()
# cat("Parameters of the component distributions:\n\t")
cat("\n")
cat("Distributions of the error components:\n")
# print(noise_params(object))
b_show(get_edist(object))
cat("\n")
invisible(object)
})
setMethod("lik_params", c(model = "MixARgen"),
function(model){
res <- callNextMethod() # get the "usual" params
## 2017-04-21 was: wrk <- mix_noiseparams(model) # get the noise params
## it seems that noise_params() is the current name of the function:
wrk <- noise_params(model) # get the noise params, TODO: option to name them?
if(length(wrk)>0)
res <- c(res, wrk)
res
})
# 2011-07-17 pravya mix_pdf etc za "MixAR" (predi byacha samo za MixARGaussian)
## cdf
setMethod("mix_cdf", signature(model="MixARgen", x="numeric", index="numeric",
xcond="missing", scale="ANY"), # note: scale is ignored here!
function(model, x, index, xcond, scale) {
## 2020-06-12 was: print("gen!")
cdf <- noise_dist(model, "cdf")
wrk <- mix_ek(model, x, index, scale=TRUE)
wrk <- cdf %of% wrk
wrk <- inner(wrk, model@prob)
wrk
})
setMethod("mix_cdf", signature(model="MixARgen", x="numeric", index="missing",
xcond="numeric", scale="ANY"),
function(model, x, index, xcond, scale) {
## 2020-06-12 was: print("gen!")
cdf <- noise_dist(model, "cdf")
wrk <- mix_hatk(model, xcond, length(xcond)+1)
# browser()
wrk <- (x - wrk)/model@scale
wrk <- cdf %of% wrk
wrk <- inner(wrk, model@prob)
wrk
})
setMethod("mix_cdf", signature(model="MixARgen", x="missing", index="missing",
xcond="numeric", scale="ANY"),
function(model, x, index, xcond, scale) {
## 2020-06-12 was: print("gen!")
cdf <- noise_dist(model, "cdf")
wrk <- mix_hatk(model, xcond, length(xcond)+1)
f <- function(x){
wrk <- (x - wrk)/model@scale
wrk <- cdf %of% wrk
wrk <- inner(wrk, model@prob)
wrk
}
f
})
## pdf
setMethod("mix_pdf", signature(model="MixARgen", x="numeric", index="numeric",
xcond="missing", scale="ANY"), # note: scale is ignored here!
function(model, x, index, xcond, scale) {
pdf <- noise_dist(model, "pdf")
wrk <- mix_ek(model, x, index, scale=TRUE)
wrk <- pdf %of% wrk
wrk <- inner(wrk, model@prob/model@scale)
wrk
})
setMethod("mix_pdf", signature(model="MixARgen", x="numeric", index="missing",
xcond="numeric", scale="ANY"),
function(model, x, index, xcond, scale) {
pdf <- noise_dist(model, "pdf")
wrk <- mix_hatk(model, xcond, length(xcond)+1)
# browser()
wrk <- (x - wrk)/model@scale
wrk <- pdf %of% wrk
wrk <- inner(wrk, model@prob/model@scale)
wrk
})
setMethod("mix_pdf", signature(model="MixARgen", x="missing", index="missing",
xcond="numeric", scale="ANY"),
function(model, x, index, xcond, scale) {
pdf <- noise_dist(model, "pdf")
wrk <- mix_hatk(model, xcond, length(xcond)+1)
# browser()
f <- function(x){
wrk <- (x - wrk)/model@scale
wrk <- pdf %of% wrk
wrk <- inner(wrk, model@prob/model@scale)
wrk
}
f
})
MixARGaussian <- setClass("MixARGaussian",
## representation(),
## prototype,
contains="MixAR"
## validity, access, where, version, sealed, package,
)
setMethod("show", "MixARGaussian",
function(object) {
callNextMethod()
cat("\n")
cat("Distributions of the error components:\n")
cat("\tstandard Gaussian\n")
cat("\n")
invisible(object)
})
## cdf
setMethod("mix_cdf", signature(model="MixARGaussian", x="numeric", index="numeric",
xcond="missing", scale="ANY"), # note: scale is ignored here!
function(model, x, index, xcond, scale) {
wrk <- mix_ek(model, x, index, scale=TRUE)
wrk <- pnorm %of% wrk
wrk <- inner(wrk, model@prob)
wrk
})
setMethod("mix_cdf", signature(model="MixARGaussian", x="numeric", index="missing", xcond="numeric",
scale="ANY"),
function(model, x, index, xcond, scale) {
wrk <- mix_hatk(model, xcond, length(xcond)+1)
# browser()
wrk <- (x - wrk)/model@scale
wrk <- pnorm %of% wrk
wrk <- inner(wrk, model@prob)
wrk
})
setMethod("mix_cdf", signature(model="MixARGaussian", x="missing", index="missing", xcond="numeric",
scale="ANY"),
function(model, x, index, xcond, scale) {
wrk <- mix_hatk(model, xcond, length(xcond)+1)
f <- function(x){
wrk <- (x - wrk)/model@scale
wrk <- pnorm %of% wrk
wrk <- inner(wrk, model@prob)
wrk
}
f
})
## pdf
setMethod("mix_pdf", signature(model="MixARGaussian", x="numeric", index="numeric", xcond="missing",
scale="ANY"), # note: scale is ignored here!
function(model, x, index, xcond, scale) {
wrk <- mix_ek(model, x, index, scale=TRUE)
wrk <- dnorm %of% wrk
wrk <- inner(wrk, model@prob/model@scale)
wrk
})
setMethod("mix_pdf", signature(model="MixARGaussian", x="numeric", index="missing", xcond="numeric",
scale="ANY"),
function(model, x, index, xcond, scale) {
wrk <- mix_hatk(model, xcond, length(xcond)+1)
# browser()
wrk <- (x - wrk)/model@scale
wrk <- dnorm %of% wrk
wrk <- inner(wrk, model@prob/model@scale)
wrk
})
setMethod("mix_pdf", signature(model="MixARGaussian", x="missing", index="missing",
xcond="numeric", scale="ANY"),
function(model, x, index, xcond, scale) {
wrk <- mix_hatk(model, xcond, length(xcond)+1)
#browser()
f <- function(x){
wrk <- (x - wrk)/model@scale
wrk <- dnorm %of% wrk
wrk <- inner(wrk, model@prob/model@scale)
wrk
}
f
})
## quantile function
setMethod("mix_qf", signature(model = "MixARGaussian", p = "numeric", x = "numeric",
index = "numeric", xcond = "missing"),
function(model, p, x, index, xcond) {
## index <- seq_along(x)[-(1:max(model@order))]
##
## TODO: verify that this is consistent with mix_pdf, mix_location, etc.
## Is this conditional quantile computed 'at time i' or 'for time i'?
## (seems it should be 'for time i', i.e. at time i-1,
## see ?mix_pdf-methods)
## ind <- max(model@order) : 1
res <- sapply(index, function(i) mix_qf(model, p, xcond = x[1:i]))
res
})
setMethod("mix_qf", signature(model = "MixARGaussian", p = "numeric", x = "missing",
index = "missing", xcond = "numeric"),
function(model, p, x, index, xcond) {
qf <- mix_qf(model, xcond = xcond)
qf(p)
})
setMethod("mix_qf", signature(model = "MixARGaussian", p = "missing", x = "missing",
index = "missing", xcond = "numeric"),
function(model, p, x, index, xcond) {
cdf <- mix_cdf(model, xcond = xcond)
f <- function(p, ..., tol = .Machine$double.eps^0.5){
if(length(p) > 1)
# sapply(p, cdf2quantile, tol = tol, ..., MoreArgs = list(cdf = cdf))
sapply(p, cdf2quantile, cdf = cdf, tol = tol, ...)
else
cdf2quantile(p, cdf, tol = tol, ...)
}
f
})
## todo: arg "1-B"
fit_mixAR <- function(x, model, init, fix, ...){
stop("There is currently no default for this function")
}
setGeneric("fit_mixAR")
setMethod("fit_mixAR", signature(x = "ANY", model = "numeric", init = "missing"),
function(x, model, init, fix, ...){
fit_mixAR(x, model, init = 1, fix = fix, ...)
})
setMethod("fit_mixAR", signature(x = "ANY", model = "numeric", init = "numeric"),
function(x, model, init, fix, ...){
model <- new("MixARGaussian", order = model)
fit_mixAR(x, model, init = init, fix = fix, ...)
})
setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "list"),
function(x, model, init, fix, ...){
wrk <- list()
for(i in init){
wrk <- c(wrk, list( fit_mixAR(x, model, init = i, fix = fix, ...) ))
}
wrk # todo: process the result, e.g. to select a best model?
})
setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "numeric"),
function(x, model, init, fix, permute, select = TRUE, ...){
if(missing(permute))
permute <- FALSE
wrk <- list()
est_shift <- missing(fix) ||
! identical(fix, "shift") # todo: tova e krapka, za da tragnat
# nestata.
templ <- est_templ(model, shift = est_shift)
for(i in 1:init){ # generate random mixAR models # todo: make 1:init safer!
ri <- em_rinit(as.numeric(x), model@order, templ)
wrk <- c(wrk, list( fit_mixAR(x, model, init = ri, fix = fix, ...) ))
# print(ri)
# print(wrk[[i]])
# browser()
}
so <- order(sapply(wrk, function(x) x$vallogf), decreasing = TRUE)
wrk <- wrk[so]
# todo: for testing, reorganise later
if(permute){ # todo: generiram otnovo vsichki permutatsii, ponezhe
# izglezhda, che inache izpuskam nesto. Obmisli i proveri pak! =>
# Otkrich kakvo stava, neveroyatno! Dvoyna greshka! permn(model@order)
# permutira elementite na order, napr. c(2, 2, 1) v rezultat
# mixAR_permute ne raboti kakto ochakvam, i ponezhe ne proveryava
# permutatsiite tichomalkom proizvezhda nevaliden model ('arcoef' ne
# saotvetsva na 'order'). Obache tova se kompensira po stechenie na
# obstoyatelstvata ot vtorata greshka, koyato se satoi v tova, che
# vmesto 'mixAR_permute' tryabva da se vika druga funktsiya (sega ste
# ya napisha). V rezultat izglezhda, che vse pak nestata se permutirat
# i rabotyat pone odnyakade kakto se ochakva kogato generiram vsichki
# permutatsii kakto po-dolu.
allperm <- permn(.nmix(model)) # todo: generate only unique ones.
# allperm <- unique( permn(model@order) ) # take only non-redundant perm.
cnt <- 0
f <- function(fm){
mo <- fm$model
locfits <- lapply(allperm,
function(p){
# locmo <- mixAR_permute(mo, p)
locmo <- mixAR_switch(mo, p)
fit_mixAR(x, model, init = locmo, fix = fix, ...)
})
cnt <<- cnt + 1
## 2020-06-12 was: cat("perm counter: ", cnt, "\n")
## print( which.max( sapply(locfits, function(z) z$vallogf) ) )
# browser() commented out on 2012-09-21
locfits[[ which.max( sapply(locfits, function(z) z$vallogf) )[1] ]]
}
newwrk <- lapply(wrk, f)
# browser()
wrk <- newwrk
}
# todo: needs more work; process the result, e.g. to select a best model;
if(isTRUE(select)){
# todo: more is needed here; the alg may stop before
# vallogf becomes NaN hence, krapka, for now.
# todo: arbitrary constant below!
loli <- sapply(wrk, function(z) if(is.nan(z$vallogf)) -Inf else z$vallogf)
ololi <- order(loli, decreasing = TRUE)
goodsigma <- sapply(wrk[ololi], function(z) all(z$model@scale > 1e-7))
if(any(goodsigma)) # todo: awful code!
res <- wrk[ololi][goodsigma][[1]]
else
res <- wrk[[ which.max(loli)[1] ]]
}else{
res <- wrk
}
res
})
## 12/09/2018 The following method now handles "raggedCoefS" objects
## note: previous content has not been changed, only added "if" conditions
setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "missing"),
function(x, model, init, fix, ...){
if(inherits(model@arcoef, "raggedCoefS")){
# if(missing(fix)) {
# mixSARfit(x, model)
# }else{
# mixSARfixfit(x, model)
# }
fix <- if(missing(fix))
FALSE
else if(is.logical(fix))
fix
## 2020-06-05 amended by Georgi
## you can't just use a completely different syntax for
## an argument that works for other methods.
## An example in inst/slowtests/example.R actually uses
## fix = "shift" and throws error (before this change).
## TODO: error message below is still confusing since the user
## may not notice that it is about mixSARfit, not mixARfit
## and even if they did, that would be equally confusing.
else if(is.character(fix) && length(fix) == 1 && fix == "shift")
TRUE
else
stop("'fix' must be boolean or missing for 'mixSARfit()'")
mixSARfit(x, model, est_shift = fix)
}else{# todo: process "..." for est_shift?
fit_mixAR(x, model, init = model, fix = fix, ...)
}# todo: do more here?
})
setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "MixAR"),
function(x, model, init, fix, ...){ # todo: process "..." for est_shift?
## 2012-11-02 premestvam processing na "fix" v mixARemFixedPoint
## todo: da machna 'fix' ot spisaka na argumentite, za da se predava s '...'.
## est_shift <- missing(fix) ||
## ! identical(fix, "shift") # todo: tova e krapka, za da tragnat
# nestata.
## todo: in this case model is a template for the result. Currently it is not
## used at all in this case but eventually should do some copying and
## other housekeeping activities. In particular, 'init' and 'model' may be
## from different subclasses of 'mixAR'.
if(missing(fix))
mixARgenemFixedPoint(x, init, ...)
else
mixARgenemFixedPoint(x, init, fix = fix, ...)
})
setMethod("fit_mixAR", signature(x = "ANY", model = "MixARGaussian", init = "MixAR"),
function(x, model, init, fix, ...){ # todo: process "..." for est_shift?
est_shift <- missing(fix) || ! identical(fix, "shift") # tova beshe krapka, za
# da tragnat nestata.
## todo: in this case model is a template for the result. Currently it is not
## used at all in this case but eventually should do some copying and
## other housekeeping activities. In particular, 'init' and 'model' may be
## from different subclasses of 'mixAR'.
mixARemFixedPoint(x, init, est_shift = est_shift, ...)
})
# build_mixAR <- function(x, ncomp = 2, pmax = 2, ...){
# # stop("There is currently no default for this funciton")
# if(ncomp < 2)
# ncomp <- 2
#
# for(g in 2:ncomp){
#
#
# }
# }
# setGeneric("build_mixAR")
#
# setMethod("build_mixAR", signature(x = "ANY", model = "numeric", init = "missing"),
# function(x, model, init, fix, ...){
# build_mixAR(x, model, init = 1, fix = fix, ...)
# })
setGeneric("mix_location", # 2012-11-08 new
function(model, x, index, xcond)
standardGeneric("mix_location"),
useAsDefault = FALSE)
setMethod("mix_location",
signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"),
function(model, x, index, xcond) {
wrk <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
wrk <- inner(wrk, model@prob)
wrk
})
setMethod("mix_location",
signature(model = "MixAR", x = "missing", index = "missing", xcond = "missing"),
function(model, x, index, xcond) {
model <- model
f <- function(xcond){
wrk <- mix_hatk(model, xcond, length(xcond)+1)
wrk <- inner(wrk, model@prob)
wrk
}
f
})
setMethod("mix_location",
signature(model = "MixAR", x = "numeric", index = "numeric", xcond = "missing"),
function(model, x, index, xcond) {
# mix_pdf and others use mix_hatk for this purpose,
# mixFilter is cleaner (and the name is better)
wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
inner(wrk, model@prob)
})
setMethod("mix_location",
signature(model = "MixAR", x = "numeric", index = "missing", xcond = "missing"),
function(model, x, index, xcond){
index <- seq_along(x)[-(1:max(model@order))]
# mix_pdf and others use mix_hatk for this purpose,
# mixFilter is cleaner (and the name is better)
wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
inner(wrk, model@prob)
})
# todo: needs to check existence of variance
setGeneric("mix_variance", # 2012-11-08 new
function(model, x, index, xcond)
standardGeneric("mix_variance"),
useAsDefault = FALSE)
setMethod("mix_variance",
signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"),
function(model, x, index, xcond) {
wrk <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
loc <- inner(wrk, model@prob)
# wrk <- inner(wrk, model@prob, function(x, y) x^2 * y )
wrk <- drop( (wrk@m^2 + model@scale^2) %*% model@prob )
wrk - loc^2
})
setMethod("mix_variance",
signature(model = "MixAR", x = "missing", index = "missing", xcond = "missing"),
function(model, x, index, xcond) {
model <- model
f <- function(xcond){
wrk <- mix_hatk(model, xcond, length(xcond)+1)
loc <- inner(wrk, model@prob)
# smenyam, ponezhe "inner" ne e dostatachno gavkava
# wrk <- inner(wrk, model@prob, function(x, y) x^2 * y)
wrk <- drop(wrk@m^2 %*% model@prob)
wrk - loc^2
}
f
})
setMethod("mix_variance",
signature(model = "MixAR", x = "numeric", index = "numeric", xcond = "missing"),
function(model, x, index, xcond) {
# mix_pdf and others use mix_hatk for this purpose,
# mixFilter is cleaner (and the name is better)
wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
loc <- inner(wrk, model@prob)
# wrk <- inner(wrk, model@prob, function(x, y) x^2 * y )
wrk <- drop(wrk@m^2 %*% model@prob)
wrk - loc^2
})
# setMethod("mix_variance",
# signature(model = "MixAR", x = "numeric", index = "missing", xcond = "missing"),
# function(model, x, index, xcond){
# index <- seq_along(x)[-(1:max(model@order))]
# # mix_pdf and others use mix_hatk for this purpose,
# # mixFilter is cleaner (and the name is better)
# wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
# inner(wrk, model@prob)
# })
setGeneric("mix_moment", # 2012-11-08 new
function(model, x, index, xcond, k)
standardGeneric("mix_moment"),
useAsDefault = FALSE)
setMethod("mix_moment",
signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"),
function(model, x, index, xcond, k) {
locall <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
loc <- inner(locall, model@prob)
gmom <- cbind(1, 0, model@scale^2, 0, 3*model@scale^4, 0) # todo: this is temp!
facbin <- choose(k, k:0)
facloc <- outer(drop(locall@m), k:0, "^")
faceps <- gmom[ , 0:k + 1]
# tuk sa obiknoveni matristi
# sum( inner(facloc * faceps * model@prob, facbin) )
sum((facloc * faceps * model@prob) %*% facbin)
})
setGeneric("mix_central_moment", # 2012-11-08 new
function(model, x, index, xcond, k)
standardGeneric("mix_central_moment"),
useAsDefault = FALSE)
setMethod("mix_central_moment", # note: k is scalar (currently)
signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"),
function(model, x, index, xcond, k) {
locall <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
loc <- inner(locall, model@prob)
# gmom <- cbind(1, 0, model@scale^2, 0, 3*model@scale^4, 0)
gmom <- outer(model@scale, 0:k, "^") * noise_moment(model, 0:k)
facbin <- choose(k, k:0)
facloc <- outer(drop(locall@m) - loc, k:0, "^")
faceps <- gmom[ , 0:k + 1]
# tuk sa obiknoveni matristi
# sum( inner(facloc * faceps * model@prob, facbin) )
sum( (facloc * faceps * model@prob) %*% facbin )
})
mix_kurtosis <- function(...){
mix_central_moment(..., k = 4) / mix_central_moment(..., k = 2)^2
}
mix_ekurtosis <- function(...){
mix_kurtosis(...) - 3
}
################
noise_moment <- function(model, k){ #elements of k are >= 0; odd moments are zero, var is 1,
if(max(k) <= 3)
c(1, 0, 1, 0)[k+1]
else
c(1, 0, 1, 0, rep(c(NA, 0), length.out = max(k-3)))[k+1]
}
setGeneric("noise_moment")
setMethod("noise_moment", signature(model = "MixARGaussian", k = "numeric"),
function(model, k){
wrk <-
if(max(k) <= 5)
c(1, 0, 1, 0, 3, 0)[k+1]
else
stdnormmoment(k)
res <- rep(wrk, .nmix(model))
if(length(k) > 1) # one row for each component
matrix(res, byrow = TRUE, nrow = .nmix(model))
else
res
})
setMethod("noise_moment", signature(model = "MixARgen", k = "numeric"),
function(model, k){
dist <- get_edist(model)
## 2021-04-30 :TODO:
## the following throw error since there is no component 'moment':
## noise_moment(exampleModels$WL_At, 3) # gives error
## noise_moment(exampleModels$WL_ibm_gen, 2) # gives error
wrk <- sapply(dist, function(x) x$moment(k))
if(length(k) > 1) # one row for each component
t(wrk) # sapply returns one column for each, so transpose
else
wrk
})
setGeneric("noise_rand",
function(model, expand = FALSE)
standardGeneric("noise_rand"),
useAsDefault = FALSE)
setMethod("noise_rand", signature(model = "MixAR"),
function(model, expand = FALSE){
mess <- paste("No method for argument of class", class(model))
stop(mess)
})
setMethod("noise_rand", signature(model = "MixARGaussian"),
function(model, expand = FALSE){
if(expand) rep(list("rnorm"), mix_ncomp(model))
else list("rnorm")
})
setMethod("noise_rand", signature(model = "MixARgen"),
function(model, expand = FALSE){
dist <- get_edist(model)
res <- lapply(dist, function(x) x[["rand"]])
if(expand) rep(res, mix_ncomp(model))
else res
})
setGeneric("noise_dist",
function(model, what, expand = FALSE)
standardGeneric("noise_dist"),
useAsDefault = FALSE)
setMethod("noise_dist", signature(model = "MixAR"),
function(model, what, expand = FALSE){
mess <- paste("No method for argument of class", class(model))
stop(mess)
})
setMethod("noise_dist", signature(model = "MixARGaussian"),
function(model, what, expand = FALSE){
res <- dist_norm[[what]]
if(expand) rep(list(res), mix_ncomp(model))
else list(res)
})
######################### new: 2011-11-16 caters for the new format of slot 'dist' in mixARgen
setGeneric("get_edist",
function(model)
standardGeneric("get_edist"),
useAsDefault = FALSE)
setMethod("get_edist", signature(model = "MixAR"),
function(model){
mess <- paste("No method for argument of class", class(model))
stop(mess)
})
setMethod("get_edist", signature(model = "MixARGaussian"),
function(model){
list(dist_norm)
})
#########################
setGeneric("noise_params",
function(model)
standardGeneric("noise_params"),
useAsDefault = FALSE)
setMethod("noise_params", signature(model = "MixAR"),
function(model){
list() # the default is no parameters (but todo: should this be numeric(0)?
})
setMethod("noise_params", signature(model = "MixARgen"), # todo: needs further thought !!
function(model){
if(!is.null(model@dist$param)){ # todo: krapka, to avoid breaking old code
## for old code
model@dist$param
}else{
dist <- get_edist(model)
sapply(dist, function(x) x$get_param())
}
})
set_noise_params <- function(model, nu){
if(!is.null(model@dist$param)){ # todo: krapka, to avoid breaking old code
model@dist$param <- nu
}else{
for(i in seq_along(model@dist))
model@dist[[i]]$set_param( nu[i] ) # very simplistic, assumes one param per dist.
}
model
}
setMethod("noise_dist", signature(model = "MixARgen"),
function(model, what, expand = FALSE){
wrk <- get_edist(model)
res <- lapply(wrk, function(x) x[[what]])
# if(is.null(res))
# stop(paste("Property", what, "has not been specified for this model."))
if(expand && length(res) == 1) rep(res, mix_ncomp(model))
else res
})
setMethod("get_edist", signature(model = "MixARgen"),
function(model){ # todo: the else part of 'if' is a patch, to avoid breaking old code
res <- if(!is.null(model@dist$generator)){
par <- model@dist$param
if(is.null(par) || length(par) == 0)
model@dist$generator()
else
model@dist$generator(par)
}else{ # tova li e za savmestimost? todo: Utochni!
model@dist
}
res
})
setGeneric("show_diff",
function(model1, model2)
standardGeneric("show_diff"),
useAsDefault=FALSE)
setMethod("show_diff", signature(model1="MixAR", model2="MixAR"),
function(model1, model2){
if(!identical(model1@order, model2@order)){
cat("The two models need to have the same number of components\n",
" and the same AR orders.")
return("")
}
g <- length(model1@prob)
p <- max(model1@order)
f <- function(i){ wrkloc <- c( model1@prob[i] - model2@prob[i]
, model1@shift[i] - model2@shift[i]
, model1@scale[i] - model2@scale[i]
, model1@order[i])
if(p > 0) c(wrkloc, model1@arcoef[[i]] - model2@arcoef[[i]])
else wrkloc
}
wrk <- lapply(1:g, f)
mcoef <- ragged2char(wrk)
colnames(mcoef) <- c("prob", "shift", "scale", "order",
if(p > 0) paste("ar_", seq_len(p), sep="")
else character(0))
rownames(mcoef) <- paste("Comp_", 1:g, sep="")
print(mcoef, na.print="", quote=FALSE, justify="right")
})
# todo: if all subclasses of mixAR can return a vector of
# descriptions will not need so many additional methods.
setMethod("show_diff", signature(model1="MixARgen", model2="MixARgen"),
function(model1, model2){
callNextMethod()
cat("\n") # todo: better formatting (in a table?)
cat("Distributions of the error components:\n")
cat("\n\tModel 1:\n"); b_show(get_edist(model1))
cat("\n\tModel 2:\n"); b_show(get_edist(model2))
invisible("")
})
setMethod("show_diff", signature(model1="MixARgen", model2="MixARGaussian"),
function(model1, model2){
callNextMethod()
cat("\n") # todo: better formatting (in a table?)
cat("Distributions of the error components:\n")
cat("\n\tModel 1\n"); b_show(get_edist(model1))
cat("\n\tModel 2\n"); cat("\t All components: Standard normal distribution\n")
invisible("")
})
setMethod("show_diff", signature(model1="MixARGaussian", model2="MixARgen"),
function(model1, model2){
callNextMethod()
cat("\n") # todo: better formatting (in a table?)
cat("Distributions of the error components:\n")
cat("\n\tModel 1\n"); cat("\t All components: Standard normal distribution\n")
cat("\n\tModel 2\n"); b_show(get_edist(model2))
invisible("")
})
# 2012-11-02 new arg. 'drop'
parameters <- function(model, namesflag = FALSE, drop = character(0)){
coef(model)
}
setGeneric("parameters")
setMethod("parameters", "MixAR",
function(model, namesflag = FALSE, drop = character(0)){ # todo: return named list!
nams <- c("order", "prob", "shift", "scale")
nams <- nams[!(nams %in% drop)] # todo: process arcoef similarly?
res <- unlist(c(lapply(nams, function(x) slot(model, x)),
model@arcoef@a
))
if(namesflag)
names(res) <- c(sapply(nams,
function(x) paste(x, 1:.nmix(model), sep="")),
paste("ar_",
rep(1:.nmix(model), times = model@order),
unlist(lapply(model@order, function(x) seq_len(x))),
sep=""))
res
})
"parameters<-" <- function(model, value){
stop("'parameters<-' has no default.")
}
setGeneric("parameters<-")
setMethod("parameters<-", "MixAR",
function(model, value){ # todo: return named list!
nams <- c("order", "prob", "shift", "scale") # todo: put this in a variable?
g <- .nmix(model)
p <- model@order # constant order
# model@order <- value[1:g] # krapka, za da raboti sas simulatsiite!
model@prob <- value[(g+1):(2*g)]
model@shift <- value[(2*g+1):(3*g)]
model@scale <- value[(3*g+1):(4*g)]
model@arcoef <- rag_modify(model@arcoef, value[-(1:(4*g))])
# browser()
model
})
# 2020-4-20
## companion_matrix and isStable now handle mixVAR objects.
companion_matrix <- function(v, ncol = length(v), nrow = ncol){
mult <- ifelse(is.matrix(v), nrow(v), 1)
if(ncol %% mult != 0) stop("Wrong dimensions for companion matrix")
if(mult > 1){
if(ncol(v) < ncol) v <- cbind(v, matrix(0, ncol = (ncol - ncol(v)), nrow = mult))
}else{
if(length(v) < ncol)
v <- c(v, rep(0, ncol - length(v)))
}
rbind(v, diag(1, nrow = nrow-mult , ncol = ncol), deparse.level = 0 ) # no labels
}
isStable <- function(x){ # todo: make generic?
cl <-inherits(x, "MixVAR")
nc <- max(x@order)
co <- if(cl) x@arcoef else x@arcoef[] # a matrix
prob <- x@prob
if(ncol(co[]) == 0) # i.i.d. mixture
return(TRUE)
f <- function(k){
if(cl){ m <- companion_matrix(co[k, ], nc * nrow(co[]))
}else m <- companion_matrix(co[k,], nc)
prob[k] * kronecker(m,m)
}
# wrk <- lapply(seq_along(x@prob), f)
wrk <- do.call(.mplus, lapply(seq_along(x@prob), f))
abs(eigen(wrk)$values[1]) < 1 # stable if maximal eigenvalue is < 1
}
mixAR_permute <- function(model, perm){ # todo: make this generic; this is only a start.
if(all(perm == seq_along(perm)))
return(model)
lapply(c("prob", "order", "shift", "scale"),
function(nam) slot(model, nam) <<- slot(model, nam)[perm] )
model@arcoef@a <- model@arcoef@a[perm]
model@arcoef@p <- model@arcoef@p[perm]
# print(model)
model
}
mixAR_switch <- function(model, perm){ # todo: check the permutation?
if(all(perm == seq_along(perm)))
return(model)
lapply(c("prob", "shift", "scale"), # note: "order" is omitted here
function(nam) slot(model, nam) <<- slot(model, nam)[perm] )
p <- model@order
m <- model@arcoef[perm]
for(i in 1:length(p)){
if(p[i]>0)
model@arcoef@a[[i]] <- m[i,1:p[i]]
}
model
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.