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
#
#
#### 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)
```

Try the distrMod package in your browser

Any scripts or data that you put into this service are public.

distrMod documentation built on Aug. 6, 2018, 9:05 a.m.