library(mixAR)
test_that("mixAR and new() for mixAR work", {
expect_output(show(new("MixARGaussian", order = c(0, 0, 0))))
new("MixARGaussian", order = c(3, 2, 1))
m2 <- new("MixARGaussian", order = c(3, 2, 1),
arcoef = matrix(c(1:3, c(1:2, 0), c(1, 0, 0)), nrow = 3, byrow = TRUE))
m3 <- new("MixARGaussian", arcoef = list(1:3, 1:2, 1))
expect_equal(m2, m3)
as(m2, "list")
## expect_identical(m2, m3) # not identical: > class(m2@arcoef@a[[1]])
## [1] "numeric"
## > class(m3@arcoef@a[[1]])
## [1] "integer"
expect_output(show(m2))
expect_output(show(exampleModels$WL_At))
expect_output(show_diff(m2, m3))
expect_output(show_diff(exampleModels$WL_At, exampleModels$WL_ibm))
expect_output(show_diff(exampleModels$WL_ibm, exampleModels$WL_At))
expect_output(show_diff(exampleModels$WL_Bt_1, exampleModels$WL_Bt_2))
expect_message(new("MixARGaussian", order = c(1, 2, 3), arcoef = list(1:3, 1:2, 1)),
"Arg's 'arcoef' and 'order' are not consistent, ignoring 'order'")
new("MixARGaussian", order = c(3, 2, 1), prob = 1)
new("MixARGaussian", order = c(3, 2, 1), prob = c(0.5, 0.25, 0.25) )
expect_error(new("MixARGaussian", order = c(3, 2, 1), prob = c(0.5, 0.25)),
"length of 'prob' should be 1 or the length of 'order'" )
new("MixARGaussian", order = c(3, 2, 1), shift = 1)
new("MixARGaussian", order = c(3, 2, 1), shift = c(0.5, 0.25, 0.25) )
expect_error(new("MixARGaussian", order = c(3, 2, 1), shift = c(0.5, 0.25)),
"length of 'shift' should be 1 or the length of 'order'" )
new("MixARGaussian", order = c(3, 2, 1), scale = 1)
new("MixARGaussian", order = c(3, 2, 1), scale = c(0.5, 0.25, 0.25) )
expect_error(new("MixARGaussian", order = c(3, 2, 1), scale = c(0.5, 0.25)),
"length of 'scale' should be 1 or the length of 'order'" )
new("MixARGaussian", model = exampleModels$WL_At)
expect_output(show_diff(new("MixARGaussian", model = exampleModels$WL_At),
exampleModels$WL_At))
expect_error(new("MixARGaussian", model = exampleModels$WL_At, order = c(2,1)),
# "Arg. 'order' (if present) must be the same as 'model@order'." )
"Arg. 'order' .if present. must be the same as 'model@order'" )
new("MixARGaussian", model = exampleModels$WL_At, prob = c(0.5, 0.5))
new("MixARGaussian", model = exampleModels$WL_At, prob = 1)
expect_error(new("MixARGaussian", model = exampleModels$WL_At, prob = c(0.5, 0.25, 0.25)),
"Arg. 'prob' .if present. must have the same length as 'model@prob'")
new("MixARGaussian", model = exampleModels$WL_At, shift = c(0.5, 0.5))
new("MixARGaussian", model = exampleModels$WL_At, shift = 1)
expect_error(new("MixARGaussian", model = exampleModels$WL_At, shift = c(0.5, 0.25, 0.25)),
"Arg. 'shift' .if present. must have the same length as 'model@shift'")
new("MixARGaussian", model = exampleModels$WL_At, scale = c(0.5, 0.5))
new("MixARGaussian", model = exampleModels$WL_At, scale = 1)
expect_error(new("MixARGaussian", model = exampleModels$WL_At, scale = c(0.5, 0.25, 0.25)),
"Arg. 'scale' .if present. must have the same length as 'model@scale'")
new("MixARGaussian", model = exampleModels$WL_At, arcoef = list(0.9, 0.4))
expect_error(new("MixARGaussian", model = exampleModels$WL_At, arcoef = list(c(0.9,0.3), 0.4)),
"Arg. 'arcoef' \\(if present\\) must be consistent with 'model@order'")
new("MixARGaussian", model = exampleModels$WL_At, arcoef = matrix(c(0.9, 0.4), ncol = 1))
new("MixARGaussian", model = exampleModels$WL_At, arcoef = matrix(c(0.9, 0.4), ncol = 1),
order = c(1, 1))
## Test function mixAR
mixAR(coef = list(prob=c(.5,.5), scale=c(1,2),
arcoef=list(.5, 1.1), shift=c(0,0), order=c(1,1)))
expect_equal(mixAR(template=c(1,1)), mixAR(coef=list(order=c(1,1))))
mixAR(m2)
mixAR(m2, list(prob = c(0.5, 0.25, 0.25)))
parameters(lm(6:10 ~ 1 + I(1:5))) # equivalent to coef() for non-mixAR objects
parameters(m2)
parameters(m2, TRUE)
v <- 1:4
expect_error(parameters(v) <- 2:5, "'parameters<-' has no default")
## exampleModels$WL_Bt_3
moT_B3 <-
new("MixARgen"
, prob = c(0.3, 0.3, 0.4)
, scale = c(2, 1, 0.5)
, shift = c(5, -5, 0)
, arcoef = list(c(0.5, 0.24), c(-0.9), c(1.5, -0.74, 0.12))
# t4, t4, t10
, dist = distlist("stdt", c(4,10), fixed = c(FALSE, TRUE), tr = c(1, 1, 2))
)
expect_output(show(moT_B3))
lik_params(moT_B3)
mix_ncomp(moT_B3)
row_lengths(moT_B3)
parameters(moT_B3)
parameters(moT_B3, TRUE)
lik_params(moT_B3)
mix_ncomp(moT_B3)
row_lengths(moT_B3)
mo <- moT_B3
val <- parameters(mo)
expect_identical(mo@prob, val[4:6]) # c(0.3, 0.3, 0.4)
val[4:6] <- c(0.5, 0.25, 0.25)
prob_old <- mo@prob
parameters(mo) <- val
expect_identical(mo@prob, val[4:6])
mo@prob <- prob_old # restore to check if anything else has changed
expect_identical(mo, moT_B3)
## demonstrate reuse of existing models
exampleModels$WL_Bt_1
## moT_C2
new("MixARgen"
, model = exampleModels$WL_Bt_1
, dist = distlist(c("stdt", "stdt", "stdnorm"), c(4,7)) # t4, t7, N(0,1)
)
## moT_C3
new("MixARGaussian", model = exampleModels$WL_Bt_1 )
expect_error(new("MixARgen" , model = exampleModels$WL_Bt_1 ,
dist = c("stdt", "stdt", "stdnorm")),
"Argument 'dist' must be a list")
new("MixARgen", model = exampleModels$WL_Bt_1)
expect_error(new("MixARgen", model = exampleModels$WL_ibm),
"Cannot set the distribution since argument 'dist' is missing")
pdf1 <- mix_pdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
## plot the predictive density
## (cdf is used to determine limits on the x-axis)
cdf1 <- mix_cdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
gbutils::plotpdf(pdf1, cdf = cdf1, lq = 0.001, uq = 0.999)
mix_pdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
mix_pdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose))
mix_cdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose))
## c(352, 357) are the last 2 values of fma::ibmclose
## 367:369 fma::ibmclose has length 369
mix_pdf(exampleModels$WL_ibm, x = c(352, 357), xcond = as.numeric(fma::ibmclose))
mix_cdf(exampleModels$WL_ibm, x = c(352, 357), xcond = as.numeric(fma::ibmclose))
mix_pdf(exampleModels$WL_ibm_gen, x = c(352, 357), xcond = as.numeric(fma::ibmclose))
mix_cdf(exampleModels$WL_ibm_gen, x = c(352, 357), xcond = as.numeric(fma::ibmclose))
mix_pdf(exampleModels$WL_ibm_gen, x = as.numeric(fma::ibmclose), index = 367:369 )
mix_cdf(exampleModels$WL_ibm_gen, x = as.numeric(fma::ibmclose), index = 367:369 )
pdf1 <- mix_pdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
cdf1 <- mix_cdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
gbutils::plotpdf(pdf1, cdf = cdf1, lq = 0.001, uq = 0.999)
qf1 <- mix_qf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
tmp <- qf1(0.05)
tmp <- qf1(1:9/10)
tmp1 <- mix_qf(exampleModels$WL_ibm, p = 1:9/10, xcond = as.numeric(fma::ibmclose))
expect_equal(tmp1, tmp)
mix_qf(exampleModels$WL_ibm, p = 0.05, x = as.numeric(fma::ibmclose), index = 367:369)
mix_qf(exampleModels$WL_ibm, p = 1:9/10, x = as.numeric(fma::ibmclose), index = 367:369)
noise_dist(exampleModels$WL_ibm_gen, "cdf")
noise_dist(exampleModels$WL_ibm_gen, "pdf")
noise_dist(exampleModels$WL_ibm_gen, "pdf", expand = TRUE)
noise_dist(exampleModels$WL_ibm_gen, "cdf", expand = TRUE)
pdf1 <- mix_pdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
cdf1 <- mix_cdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose))
gbutils::plotpdf(pdf1, cdf = cdf1, lq = 0.001, uq = 0.999)
pdf1gen <- mix_pdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose))
cdf1gen <- mix_cdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose))
gbutils::plotpdf(pdf1gen, cdf = cdf1gen, lq = 0.001, uq = 0.999)
length(fma::ibmclose)
cdf1gena <- mix_cdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose)[-(369:369)])
pdf1gena <- mix_pdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose)[-(369:369)])
gbutils::plotpdf(pdf1gena, cdf = cdf1gena, lq = 0.001, uq = 0.999)
pdf1a <- mix_pdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose)[-(369:369)])
cdf1a <- mix_cdf(exampleModels$WL_ibm, xcond = as.numeric(fma::ibmclose)[-(369:369)])
gbutils::plotpdf(pdf1a, cdf = cdf1a, lq = 0.001, uq = 0.999)
cdf1gena <- mix_cdf(exampleModels$WL_ibm_gen, xcond = as.numeric(fma::ibmclose)[-(369:369)])
# cond_loglik(exampleModels$WL_ibm, fma::ibmclose) # doesn't work currently, no method for `ts'
cond_loglik(exampleModels$WL_ibm, as.numeric(fma::ibmclose))
cond_loglik(exampleModels$WL_ibm_gen, as.numeric(fma::ibmclose))
ts1gen <- mixAR_sim(exampleModels$WL_ibm_gen, n = 30, init = c(346, 352, 357), nskip = 0)
plot(ts1gen)
plot(mixAR_sim(exampleModels$WL_ibm_gen, n = 100, init = c(346, 352, 357), nskip = 0),
type = "l")
plot(diff(mixAR_sim(exampleModels$WL_ibm_gen, n = 100, init = c(346, 352, 357), nskip = 0)),
type = "l")
noise_dist(exampleModels$WL_ibm_gen, "Fscore")
prob <- exampleModels$WL_ibm@prob
scale <- exampleModels$WL_ibm@scale
arcoef <- exampleModels$WL_ibm@arcoef@a
mo_WLt3 <- new("MixARgen", prob = prob, scale = scale, arcoef = arcoef,
dist = list(fdist_stdt(3)))
mo_WLt30 <- new("MixARgen", prob = prob, scale = scale, arcoef = arcoef,
dist = list(fdist_stdt(30)))
f <- make_fcond_lik(exampleModels$WL_ibm, as.numeric(fma::ibmclose))
f(as.numeric(fma::ibmclose))
f2 <- make_fcond_lik(exampleModels$WL_At, as.numeric(fma::ibmclose))
f2(as.numeric(fma::ibmclose))
## data(ibmclose, package = "fma") # `ibmclose'
ibmclose <- as.numeric(fma::ibmclose)
length(ibmclose) # 369
max(exampleModels$WL_ibm@order) # 2
## compute point predictions for t = 3,...,369
mix_location(exampleModels$WL_ibm, ibmclose)
## compute one-step point predictions for t = 360,...369
mix_location(exampleModels$WL_ibm, ibmclose, index = 369 - 9:0 )
f <- mix_location(exampleModels$WL_ibm) # a function
## predict the value after the last
f(ibmclose)
## a different way to compute one-step point predictions for t = 360,...369
sapply(369 - 10:1, function(k) f(ibmclose[1:k]))
## the results are the same, but notice that xcond gives past values
## while index above specifies the times for which to compute the predictions.
identical(sapply(369 - 10:1, function(k) f(ibmclose[1:k])),
mix_location(exampleModels$WL_ibm, ibmclose, index = 369 - 9:0 ))
## conditional variance
f <- mix_variance(exampleModels$WL_ibm) # a function
## predict the value after the last
f(ibmclose)
## a different way to compute one-step point predictions for t = 360,...369
sapply(369 - 10:1, function(k) f(ibmclose[1:k]))
## the results are the same, but notice that xcond gives past values
## while index above specifies the times for which to compute the predictions.
identical(sapply(369 - 10:1, function(k) f(ibmclose[1:k])),
mix_variance(exampleModels$WL_ibm, ibmclose, index = 369 - 9:0 ))
# interesting example
# bimodal distribution, low kurtosis, 4th moment not much larger than 2nd
moWL <- exampleModels$WL_ibm
mix_location(moWL,xcond = c(500,450))
mix_kurtosis(moWL,xcond = c(500,450))
mix_ekurtosis(moWL,xcond = c(500,450))
f1pdf <- mix_pdf(moWL,xcond = c(500,450))
f1cdf <- mix_cdf(moWL,xcond = c(500,450))
gbutils::plotpdf(f1pdf,cdf=f1cdf)
gbutils::plotpdf(f1cdf,cdf=f1cdf)
f1cdf(c(400,480))
mix_variance(moWL,xcond = c(500,450))
mix_central_moment(moWL,xcond = c(500,450), k=2)
mix_moment(moWL,xcond = c(500,450), k=2)
sqrt(mix_variance(moWL,xcond = c(500,450)))
sqrt(mix_central_moment(moWL,xcond = c(500,450), k=2))
bycomp <- list(list(0.1, 10, 0.11, 1),
list(0.2, 20, c(0.11, 0.22), 2),
list(0.3, 30, c(0.11, 0.22, 0.33), 3) )
bytype <- tomarparambyType(bycomp)
identical(bycomp, tomarparambyComp(bytype)) # TRUE
permuteArpar(bycomp)
companion_matrix(4:1)
companion_matrix(4:1, ncol=6)
companion_matrix(4:1, ncol=6, nrow=3)
isStable(exampleModels$WL_I)
isStable(exampleModels$WL_II)
stdt3v <- fdist_stdt(3, fixed = FALSE)
companion_matrix(11:14)
companion_matrix(11:15, 6)
expect_warning(companion_matrix(11:15, ncol = 2)) # ncol < length(v)
isStable(exampleModels$WL_I)
isStable(exampleModels$WL_II)
isStable(new("MixARGaussian", order = c(0, 0, 0)))
## missing components are filled with 'filler', extended accordingly
mixAR:::.canonic_coef(list(order = c(2,3)), filler = NA)
# here 'scale' is replicated, the missing 'shift' is inserted
mo <- list(order = c(2,3), prob = c(0.4, 0.6), scale = 1,
arcoef = list(c(0.5, -0.5), c(1.1, 0.0, -0.5)) )
mixAR:::.canonic_coef(mo, filler = NA)
### seasonal models
probS <- c(0.5, 0.5)
sigmaS <- c(5, 1)
ar1 <- list(c(0.5, -0.5), c(1.1, 0, -0.5))
ar12 <- list(0, c(-0.3, 0.1))
s <- 12
rag <- new("raggedCoefS", a=ar1, as=ar12, s=s)
modelS <- new("MixARGaussian", prob=probS, scale=sigmaS, arcoef=rag)
yS <- mixAR_sim(modelS, n=500, init=rep(0,24))
expect_error(new("MixARGaussian", prob=probS, scale=sigmaS,
arcoef=list(a=ar1, s=s)))
new("MixARGaussian", prob=probS, scale=sigmaS, arcoef=list(ar1, ar12, s=s))
new("MixARGaussian", prob=probS, scale=sigmaS, arcoef=list(ar1, as=ar12, s=s))
new("MixARGaussian", prob=probS, scale=sigmaS, arcoef=list(a=ar1, ar12, s=s))
new("MixARGaussian", prob=probS, scale=sigmaS,
arcoef=list(a=rbind(c(0.5, -0.5, 0), c(1.1, 0, -0.5)),
as=rbind(c(0,0), c(-0.3, 0.1)) , s=s))
expect_error(new("MixARGaussian", prob=probS, scale=sigmaS,
arcoef=list(a=rbind(c(0.5, -0.5, 0), c(1.1, 0, -0.5)),
as=matrix(c(0,0), nrow=1) , s=s)))
expect_output(show(modelS))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.