inst/scripts/censoredPois.R

In distrMod: Object Oriented Implementation of Probability Models

##########################################################
### a new model -- a censored Poisson distribution
### i.e. we only observe values if they are larger than
##  a lower truncation point
##########################################################
require(distrMod)
options("newDevice"=TRUE)

CensoredPoisFamily <- function(lambda = 1, trunc.pt = 2, trafo=1){
    ## name
    name <- "Censored Poisson family"
    ## central distribution
    distribution <- Truncate(Pois(lambda = lambda), lower= trunc.pt )
    param0 <- lambda
    names(param0) <- "lambda"
    ## parameter definition
    param <- ParamFamParameter(name = "positive mean",
                               main = param0)

    ## mapping theta -> P_theta
    modifyParam <- function(theta){
                      P <- Pois(lambda = theta)
                      Truncate(P, lower = trunc.pt)}

    ## search interval for reasonable parameters
    startPar <- function(x,...) c(.Machine$double.eps,max(x))

    ## what to do in case of leaving the parameter domain
    makeOKPar <- function(param) max(param,.Machine$double.eps^.4)

    ## mapping theta -> Lambda_theta
    L2deriv.fct <- function(param){
                   lambda <- main(param)
                   fct <- function(x){}
                   body(fct) <- substitute({
                                 x/lambda-ppois(trunc.pt-1, lambda = lambda,
                                                lower.tail=FALSE)/
                                          ppois(trunc.pt, lambda = lambda,
                                                lower.tail=FALSE)},
                                list(lambda = lambda))
                   return(fct)}

    res <- L2ParamFamily(name = name, distribution = distribution,
                         param = param, modifyParam = modifyParam,
                         L2deriv.fct = L2deriv.fct,
                         startPar = startPar, makeOKPar = makeOKPar)

    ## a simplified call
    res@fam.call <- substitute(CensoredPoisFamily(lambda = l, trunc.pt = t),
                               list(l = lambda, t = trunc.pt))
    return(res)
}

## assign the model
CP <- CensoredPoisFamily(3,2)

## some observations
CP.data <- r(CP)(40)

## MLE
(m<- MLEstimator(CP.data, CP))

#Evaluations of Maximum likelihood estimate:
#-------------------------------------------
#An object of class “Estimate”
#generated by call
#  MLEstimator(x = CP.data, ParamFamily = CP)
#samplesize:   40
#estimate:
#    lambda
#  2.6221186
# (0.2737858)
#asymptotic (co)variance (multiplied with samplesize):
#[1] 2.998346
#Criterion:
#negative log-likelihood
#               59.45028

confint(m)

#A[n] asymptotic (CLT-based) confidence interval:
#          2.5 %   97.5 %
#lambda 2.085508 3.158729
#Type of estimator: Maximum likelihood estimate
#samplesize:   40
#Call by which estimate was produced:
#MLEstimator(x = CP.data, ParamFamily = CP)

plot(profile(m))

## MDE
(md.kolm<- MDEstimator(CP.data, CP))

#Evaluations of Minimum Kolmogorov distance estimate:
#----------------------------------------------------
#An object of class “Estimate”
#generated by call
#  MDEstimator(x = CP.data, ParamFamily = CP)
#samplesize:   40
#estimate:
#  lambda
#2.756499
#Criterion:
#Kolmogorov distance
#         0.04190892


(md.CvM<-  MDEstimator(CP.data, CP, distance = CvMDist,
           asvar.fct = distrMod:::.CvMMDCovariance))

#Evaluations of Minimum CvM distance estimate:
#---------------------------------------------
#An object of class “Estimate”
#generated by call
#  MDEstimator(x = CP.data, ParamFamily = CP, distance = CvMDist,
#    asvar.fct = distrMod:::.CvMMDCovariance)
#samplesize:   40
#estimate:
#    lambda
#  2.701491
# (1.926043)
#asymptotic (co)variance (multiplied with samplesize):
#[1] 148.3857
#Criterion:
#CvM distance
#  0.03700009


confint(md.CvM)

#A[n] asymptotic (CLT-based) confidence interval:
#           2.5 %   97.5 %
#lambda -1.073484 6.476466
#Type of estimator: Minimum CvM distance estimate
#samplesize:   40
#Call by which estimate was produced:
#MDEstimator(x = CP.data, ParamFamily = CP, distance = CvMDist,
#    asvar.fct = distrMod:::.CvMMDCovariance)

###--->PROBLEM:
plot(profile(md.CvM))

if(require(ROptEst)){
CP.data0 <- r(CP)(10000)
CP.data1 <- CP.data0; CP.data1[sample(1:100,10)] <- NA
print(md.ropt<- roptest(CP.data0, CP, eps=0.1, initial.est=md.CvM))
confint(md.ropt,symmetricBias())
}

#Evaluations of 1-step estimate:
#-------------------------------
#An object of class “Estimate”
#generated by call
#  roptest(x = CP.data0, L2Fam = CP, eps = 0.1, initial.est = md.CvM)
#samplesize:   10000
#estimate:
#     lambda
#  2.87822794
# (0.01379214)
#asymptotic (co)variance (multiplied with samplesize):
#[1] 1.902231
#Infos:
#     method    message
#[1,] "roptest" "1-step estimate for Censored Poisson family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 16.28873
#(partial) influence curve:
#An object of class “ContIC”
#### name:        IC of contamination type
#
#### L2-differentiable parametric family:         Censored Poisson family
#### param:      An object of class "ParamFamParameter"
#name:   positive mean
#lambda: 2.70149127445287
#
#### neighborhood radius:         10
#
#### clip:       [1] 1.628873
#### cent:       [1] -99.29597
#### stand:
#         [,1]
#[1,] 267.2244
#
#### Infos:
#     method  message
#[1,] "optIC" "optimally robust IC for ‘asMSE’"
#steps:
#[1] 1
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
#          2.5 %   97.5 %
#lambda 2.626539 3.129917
#Type of estimator: 1-step estimate
#samplesize:   10000
#Call by which estimate was produced:
#roptest(x = CP.data0, L2Fam = CP, eps = 0.1, initial.est = md.CvM)