# NOT RUN {
data(accidents)
# convert the data to vector:
accidents.obs <- unlist(accidents)
# generate MLE function:
MLE.pois <- function(dat) mean(dat)
# generate function needed for estimating the j^th moment of the
# mixing distribution via Hankel.method "explicit"
explicit.pois <- function(dat, j){
mat <- matrix(dat, nrow = length(dat), ncol = j) -
matrix(0:(j-1), nrow = length(dat), ncol = j, byrow = TRUE)
return(mean(apply(mat, 1, prod)))
}
# construct a 'datMix' object:
accidents.dM <- datMix(accidents.obs, dist = "pois", discrete = TRUE,
Hankel.method = "explicit",
Hankel.function = explicit.pois,
theta.bound.list = list(lambda = c(0, Inf)),
MLE.function = MLE.pois)
# define the penalty:
pen <- function(j, n) j * log(n)
## complexity estimation:
set.seed(0)
res <- paramHankel(accidents.dM, j.max = 5, B = 1000, ql = 0.025, qu = 0.975)
# plot the results:
plot(res, breaks = 8, ylim = c(0, 0.8))
# }
Run the code above in your browser using DataLab