## ---- include = FALSE--------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----setup-------------------------------------------------------------------- library(spINAR) ## ----------------------------------------------------------------------------- # Load the package library(spINAR) ## ----------------------------------------------------------------------------- ice <- c(0, 0, 0, 1, 1, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 2, 3, 2, 1, 0, 1, 1, 0, 2, 0, 0, 1, 2, 2, 1, 1, 1, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 2, 1, 1, 2, 2, 2, 1, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 0, 5, 1, 2, 0, 0, 0, 1, 2, 0, 1, 1, 3, 2, 3, 1, 2, 3, 2, 1, 1, 1, 1, 2, 0, 1, 0, 0, 1, 1, 1, 2, 3, 1, 1, 1, 2, 2, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 0, 1, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 2, 1, 0, 1, 0, 0, 1, 0, 2, 1, 1, 2, 0, 1, 0, 0, 2, 0, 1, 1, 0, 1, 1, 0, 1, 0, 3, 1, 1, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 2, 2, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 2, 2, 2, 2, 2, 0, 2, 2, 0, 2, 2, 2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 3, 4, 4, 4, 3, 3, 3, 2, 1, 2, 2, 2, 2, 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 0, 1, 2, 2, 1, 1, 1, 0, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2, 1, 1, 0, 1, 1, 2, 3, 2, 2, 2, 1, 2, 1, 0, 1, 1, 1, 0, 0, 4, 3, 2, 1, 3, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 1, 1, 2, 1, 0, 0, 0, 0, 0, 1, 2, 2, 2, 0, 0, 0, 1, 3, 1, 2, 1, 0, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 2, 0, 1, 1, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 1, 0, 0, 1, 1, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 1, 1, 1, 1, 1, 0, 1, 1, 2, 2, 2, 1, 1, 0, 1, 0, 1, 0, 1, 2, 3, 2, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 1, 2, 2, 0, 1, 0, 0, 1, 0, 1, 2, 1, 1, 1, 0, 0, 2, 0, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 2, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 3, 3, 5, 1, 2, 2, 2, 2, 3, 2, 3, 3, 2, 3, 3, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 3, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 3, 3, 3, 2, 1, 2, 1, 0, 1, 2, 1, 2, 1, 1, 1, 2, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 2, 2, 4, 3, 2, 0, 2, 2, 2, 0, 2, 2, 2, 2, 0, 1, 2, 4, 2, 1, 0, 1, 3, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 2, 1, 1, 2, 2, 2, 0, 2, 1, 2, 0, 2, 2, 0, 0, 0, 2, 2, 1, 2, 2, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 2, 3, 2, 2, 0, 2, 1, 1, 2, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 1, 1, 1, 0, 2, 1, 1, 1, 1, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 2, 1, 2, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 3, 2, 2, 0, 0, 1, 1, 0, 2, 0, 1, 0, 1, 1, 2, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 2, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 4, 1, 2, 1, 1, 3, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 4, 0, 2, 2, 3, 1, 1, 1, 1, 2, 0, 3, 3, 2, 1, 1, 1, 1, 0, 2, 2, 2, 1, 0, 1, 0, 0, 1, 2, 1, 2, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 2, 0, 0, 2, 0, 4, 2, 2, 0, 1, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 3, 4, 5, 4, 4, 3, 2, 2, 2, 2, 2, 1, 1, 2, 1, 2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 2, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 2, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 1, 0, 3, 1, 2, 3, 3, 1, 1, 0, 0, 2, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 0, 0, 0, 1, 1, 1, 1, 1, 1, 3, 1, 1, 2, 0, 0, 2, 1, 1, 2, 3, 2, 7, 3, 6, 5, 0, 1, 1, 1, 1, 2, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 2, 0, 1, 1, 1, 1, 1, 1, 2, 0, 1, 2, 1, 1, 2, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2, 1, 2, 3, 1, 1, 1, 0, 1, 1, 1, 2, 1, 2, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 2, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) ## ----------------------------------------------------------------------------- # p=1 xmax <- max(ice) est_sp_1 <- spinar_est(ice,1) mu_e_sp_1 <- sum(est_sp_1[-1]*(0:xmax)) sigma2_e_sp_1 <- sum(est_sp_1[-1]*(0:xmax)^2) - mu_e_sp_1^2 c(est_sp_1, mu_e_sp_1,sigma2_e_sp_1) # p=2 est_sp_2 <- spinar_est(ice,2) mu_e_sp_2 <- sum(est_sp_2[-(1:2)]*(0:xmax)) sigma2_e_sp_2 <- sum(est_sp_2[-(1:2)]*(0:xmax)^2) - mu_e_sp_2^2 c(est_sp_1, mu_e_sp_2,sigma2_e_sp_2) ## ----------------------------------------------------------------------------- # p=1 est_p_1 <- spinar_est_param(ice,1,"mom","poi") est_p_1_v <- c(est_p_1[1],dpois(0:xmax,est_p_1[2])) unname(c(est_p_1_v, est_p_1[2], est_p_1[2])) # p=2 est_p_2 <- spinar_est_param(ice,2,"mom","poi") est_p_2_v <- c(est_p_2[1:2],dpois(0:xmax,est_p_2[3])) unname(c(est_p_2_v, est_p_2[3], est_p_2[3])) ## ----------------------------------------------------------------------------- # Setting of parametrizations (Jentsch and Weiß, 2017) B <- 10^4 level <- 0.05 M <- 10 mu_est <- mean(ice) ## ----------------------------------------------------------------------------- # p=1 xstar_p <- spinar_boot(ice, 1, B, "p", "mom", "poi", progress = FALSE)$x_star mu_est_star <- apply(xstar_p, 2, mean) mu_est_star_cent <- mu_est_star - mu_est srt <- sort(mu_est_star_cent) c(mu_est - quantile(srt,1-level,names=FALSE), mu_est - quantile(srt,level,names=FALSE)) # p=2 xstar_p <- spinar_boot(ice, 2, B, "p", "mom", "poi", progress = FALSE)$x_star mu_est_star <- apply(xstar_p, 2, mean) mu_est_star_cent <- mu_est_star - mu_est srt <- sort(mu_est_star_cent) c(mu_est - quantile(srt,1-level,names=FALSE), mu_est - quantile(srt,level,names=FALSE)) ## ----eval=FALSE--------------------------------------------------------------- # # p=1 # mu_est_cent <- mu_e_sp_1/(1-est_sp_1[1]) # xstar_sp <- spinar_boot(ice, 1, B, "sp", progress = FALSE)$x_star # mu_est_star <- apply(xstar_sp, 2, mean) # mu_est_star_cent <- mu_est_star - mu_est_cent # srt <- sort(mu_est_star_cent) # c(mu_est - quantile(srt,1-level,names=FALSE), mu_est - quantile(srt,level,names=FALSE)) # # # p=2 # mu_est_cent <- mu_e_sp_2/(1-est_sp_2[1]-est_sp_2[2]) # xstar_sp <- spinar_boot(ice, 2, B, "sp", progress = FALSE)$x_star # mu_est_star <- apply(xstar_sp, 2, mean) # mu_est_star_cent <- mu_est_star - mu_est_cent # srt <- sort(mu_est_star_cent) # c(mu_est - quantile(srt,1-level,names=FALSE), mu_est - quantile(srt,level,names=FALSE)) ## ----------------------------------------------------------------------------- # Load the data set data <- c(1,1,0,2,1,4,4,5,4,0,2,1,0,0,1,2,1,3,1,0,1,2,1,0,0,1,0,1,0,0,2,0,2,2,0, 2,1,0,1,2,1,1,0,0,0,0,0,0,1,2,2) ## ----------------------------------------------------------------------------- # Unpenalized estimation of the innovation distribution est_unpenal <- spinar_est(data,1) # Penalized estimation of the innovation distribution est_penal <- spinar_penal(data,1,0,1) ## ----eval=FALSE--------------------------------------------------------------- # par(mfrow=c(1,2)) # barplot(est_unpenal[-1], ylim=c(0,1),names.arg=0:5, # main="Unpenalized estimated \n innovation distribution") # barplot(est_penal[-1], ylim=c(0,1),names.arg=0:5, # main="Penalized estimated \n innovation distribution")