Nothing
##########################################################
### 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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.