Nothing
inst/scripts/BinomialModel.R
In ROptEst: Optimally Robust Estimation
###############################################################################
## Example: Binomial Family
###############################################################################
require(ROptEst)
options("newDevice"=TRUE)
#-------------------------------------------------------------------------------
## Preparations
#-------------------------------------------------------------------------------
## generates Binomial Family with
## m = 25 and probability of success theta = 0.25
B <- BinomFamily(size = 25, prob = 0.25)
B # show B
#An object of class "BinomFamily"
#### name: Binomial family
#
#### distribution: Distribution Object of Class: Binom
# size: 25
# prob: 0.25
#
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### props:
#[1] "The Binomial family is symmetric with respect to prob = 0.5;"
#[2] "i.e., d(Binom(size, prob))(k)=d(Binom(size,1-prob))(size-k)"
plot(B) # plot of Binom(size = 25, prob = 0.25) and L_2 derivative
checkL2deriv(B)
#precision of centering: -3.103725e-15
#precision of Fisher information:
# prob
#prob -2.842171e-14
#$maximum.deviation
#[1] 2.842171e-14
#-------------------------------------------------------------------------------
## classical optimal IC
#-------------------------------------------------------------------------------
IC0 <- optIC(model = B, risk = asCov())
IC0 # show IC
#An object of class IC
#### name: Classical optimal influence curve for Binomial family
#### L2-differentiable parametric family: Binomial family
#
#### 'Curve': An object of class EuclRandVarList
#Domain: Real Space with dimension 1
#[[1]]
#length of Map: 1
#Range: Real Space with dimension 1
#
#### Infos:
# method message
#[1,] "optIC" "optimal IC in sense of Cramer-Rao bound"
plot(IC0) # plot IC
checkIC(IC0)
#precision of centering: -2.403827e-17
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
Risks(IC0)
#$asCov
# prob
#prob 0.0075
#
#$trAsCov
#[1] 0.0075
#-------------------------------------------------------------------------------
## lower case radius
#-------------------------------------------------------------------------------
lowerCaseRadius(L2Fam = B, neighbor = ContNeighborhood(), risk = asMSE())
#lower case radius
# 1.197809
lowerCaseRadius(L2Fam = B, neighbor = TotalVarNeighborhood(), risk = asMSE())
#lower case radius
# 1.130615
#-------------------------------------------------------------------------------
## L_2 family + infinitesimal neighborhood
#-------------------------------------------------------------------------------
RobB1 <- InfRobModel(center = B, neighbor = ContNeighborhood(radius = 0.5))
RobB1 # show RobB1
#An object of class InfRobModel
####### center: An object of class "BinomFamily"
#### name: Binomial family
#
#### distribution: Distribution Object of Class: Binom
# size: 25
# prob: 0.25
#
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### props:
#[1] "The Binomial family is symmetric with respect to prob = 0.5;"
#[2] "i.e., d(Binom(size, prob))(k)=d(Binom(size,1-prob))(size-k)"
#
####### neighborhood: An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
(RobB2 <- InfRobModel(center = B, neighbor = TotalVarNeighborhood(radius = 0.5)))
#An object of class InfRobModel
####### center: An object of class "BinomFamily"
#### name: Binomial family
#
#### distribution: Distribution Object of Class: Binom
# size: 25
# prob: 0.25
#
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### props:
#[1] "The Binomial family is symmetric with respect to prob = 0.5;"
#[2] "i.e., d(Binom(size, prob))(k)=d(Binom(size,1-prob))(size-k)"
#
####### neighborhood: An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
#-------------------------------------------------------------------------------
## OBRE solution
#-------------------------------------------------------------------------------
system.time(ICA <- optIC(model=RobB1, risk=asAnscombe(),
verbose=TRUE,lower=NULL,upper=10))
#-------------------------------------------------------------------------------
## MSE solution
#-------------------------------------------------------------------------------
system.time(IC1 <- optIC(model=RobB1, risk=asMSE()))
# user system elapsed
# 3.62 0.00 3.62
IC1
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1213661
#### cent: [1] -0.003061948
#### stand:
# prob
#prob 0.01222674
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
checkIC(IC1)
#precision of centering: -3.17283e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
Risks(IC1)
#$asCov
# prob
#prob 0.008544305
#
#$asBias
#$asBias$value
#[1] 0.1213661
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008544305
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01222674
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
getRiskIC(IC1, asBias(), ContNeighborhood()) # standardized bias
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) convex contamination neighborhood"
#
#$asBias$value
#[1] 0.1213661
getRiskIC(IC1, asMSE(), ContNeighborhood(radius = 0.5))
#$asMSE
#$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$neighborhood
#[1] "(uncond.) convex contamination neighborhood with radius 0.5"
#
#$asMSE$radius
#[1] 0.5
(Cov1 <- getRiskIC(IC1, asCov()))
#$asCov
#$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$value
#$asCov$value[[1]]
#[1] 0.008544305
#
#$asCov$value$asCov
#$asCov$value$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$value$asCov$value
# prob
#prob 0.008544305
(mse1 <- getRiskIC(IC1, asMSE(), TotalVarNeighborhood(radius = 0.5)))
#$asMSE
#$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$radius
#[1] 0.5
#
#$asMSE$value
#[1] 0.01222674
(bias1 <- getRiskIC(IC1, asBias(), TotalVarNeighborhood()))
#
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$value
#[1] 0.1213661
## only suboptimal -> ToDo-List
addRisk(IC1) <- list(Cov1, mse1, bias1)
Risks(IC1)
#$asCov
#$asCov[[1]]
#[1] 0.008544305
#
#$asCov$asCov
#$asCov$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value
# prob
#prob 0.008544305
#
#
#$asCov$asCov
#$asCov$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value
#$asCov$asCov$value[[1]]
#[1] 0.008544305
#
#$asCov$asCov$value$asCov
#$asCov$asCov$value$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value$asCov$value
# prob
#prob 0.008544305
#
#
#$asBias
#$asBias$value
#[1] 0.1213661
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#$asBias$asBias
#$asBias$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$asBias$value
#[1] 0.1213661
#
#
#$asBias$asBias
#$asBias$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$asBias$value
#[1] 0.1213661
#
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008544305
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01222674
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
#
#$asMSE$asMSE
#$asMSE$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$asMSE$radius
#[1] 0.5
#
#$asMSE$asMSE$value
#[1] 0.01222674
#
#
#$asMSE$asMSE
#$asMSE$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$asMSE$radius
#[1] 0.5
#
#$asMSE$asMSE$value
#[1] 0.01222674
plot(IC1)
system.time(IC2 <- optIC(model=RobB2, risk=asMSE()))
# user system elapsed
# 9.46 0.02 9.47
IC2
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clipLo: [1] -0.1081481
#### clipUp: [1] 0.1158377
#### stand:
# prob
#prob 0.02208611
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
checkIC(IC2)
#precision of centering: -1.745301e-13
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.745301e-13
Risks(IC2)
#$asCov
# prob
#prob 0.03342380
#
#$asBias
#$asBias$value
#[1] 0.2239858
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.03342380
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.04596621
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
#
getRiskIC(IC2, asMSE(), TotalVarNeighborhood(radius = 0.5))
#$asMSE
#$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$radius
#[1] 0.5
#
#$asMSE$value
#[1] 0.04596621
getRiskIC(IC2, asBias(), TotalVarNeighborhood())
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$value
#[1] 0.2239858
getRiskIC(IC2, asBias(), ContNeighborhood())
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) convex contamination neighborhood"
#
#$asBias$value
#[1] 0.2239858
Cov2 <- getRiskIC(IC2, asCov())
addRisk(IC2) <- Cov2
Risks(IC2)
#$asCov
#$asCov[[1]]
#[1] 0.03342380
#
#$asCov$asCov
#$asCov$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value
# prob
#prob 0.03342380
#
#$asBias
#$asBias$value
#[1] 0.2239858
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.03342380
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.04596621
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
plot(IC2)
#-------------------------------------------------------------------------------
## lower case solutions
#-------------------------------------------------------------------------------
(IC3 <- optIC(model=RobB1, risk=asBias()))
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1098910
#### cent: [1] -1.333333
#### stand:
# [,1]
#[1,] 1
#### lowerCase: [1] -0.3316026
#
#### Infos:
# method message
#[1,] "optIC" "minimum asymptotic bias (lower case) solution"
checkIC(IC3)
#precision of centering: -1.714470e-17
#precision of Fisher consistency:
# prob
#prob 2.220446e-16
#maximum deviation
# 2.220446e-16
Risks(IC3)
#$asBias
#$asBias$value
#[1] 0.1098910
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$asCov
#[1] 0.01011106
#
#$trAsCov
#$trAsCov$value
#[1] 0.01011106
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
#[1] 0.01086581
#
#$asMSE$r
#[1] 0.25
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
plot(IC3)
(IC4 <- optIC(model=RobB2, risk=asBias()))
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
### param: An object of class "ParamFamParameter"
##name: probability of success
#prob: 0.25
##fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clipLo: [1] -0.094766
#### clipUp: [1] 0.1211501
#### stand:
# [,1]
#[1,] 1
#
#### Infos:
# method message
#[1,] "optIC" "minimum asymptotic bias (lower case) solution"
checkIC(IC4)
#precision of centering: -1.076010e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
Risks(IC4)
#$asBias
#$asBias$value
#[1] 0.2159161
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$asCov
#[1] 0.01148091
#
#$trAsCov
#$trAsCov$value
#[1] 0.01148091
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
#[1] 0.01439465
#
#$asMSE$r
#[1] 0.25
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
plot(IC4)
#-------------------------------------------------------------------------------
## Hampel solution
#-------------------------------------------------------------------------------
(IC5 <- optIC(model=RobB1, risk=asHampel(bound=clip(IC1))))
#minimal bound: 0.1098910
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1213661
#### cent: [1] -0.003061961
#### stand:
# prob
#prob 0.01222672
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for 'asHampel' with bound = 0.121"
checkIC(IC5)
#precision of centering: -3.524202e-17
#precision of Fisher consistency:
# prob
#prob -3.56339e-07
#maximum deviation
# 3.56339e-07
Risks(IC5)
#$asCov
# prob
#prob 0.008544296
#
#$asBias
#$asBias$value
#[1] 0.1213661
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008544296
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01222673
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
plot(IC5)
(IC6 <- optIC(model=RobB2, risk=asHampel(bound=Risks(IC2)$asBias$value), maxiter = 200))
#minimal bound: 0.2159161
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clipLo: [1] -0.1081481
#### clipUp: [1] 0.1158378
#### stand:
# prob
#prob 0.02208608
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for 'asHampel' with bound = 0.224"
checkIC(IC6)
#precision of centering: -1.745225e-13
#precision of Fisher consistency:
# prob
#prob -1.166564e-07
#maximum deviation
# 1.166564e-07
Risks(IC6)
#$asCov
# prob
#prob 0.03342372
#
#$asBias
#$asBias$value
#[1] 0.2239858
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.03342372
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.04596613
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
plot(IC6)
#-------------------------------------------------------------------------------
## radius minimax IC
#-------------------------------------------------------------------------------
system.time(IC7 <- radiusMinimaxIC(L2Fam=B, neighbor=ContNeighborhood(),
risk=asMSE(), loRad=0, upRad=1))
# user system elapsed
# 39.39 0.02 39.45
IC7
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3913008
#
#### clip: [1] 0.1269559
#### cent: [1] -0.004305063
#### stand:
# prob
#prob 0.01074017
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 1]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.391"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.103"
checkIC(IC7)
#precision of centering: -3.056847e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
Risks(IC7)
#$asCov
# prob
#prob 0.008272273
#
#$asBias
#$asBias$value
#[1] 0.1269559
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008272273
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01074017
#
#$asMSE$r
#[1] 0.3913008
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.3913008
plot(IC7)
system.time(IC8 <- radiusMinimaxIC(L2Fam=B, neighbor=TotalVarNeighborhood(),
risk=asMSE(), loRad=0, upRad=1))
# user system elapsed
# 163.20 0.12 168.21
IC8
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.2661058
#
#### clipLo: [1] -0.1157700
#### clipUp: [1] 0.1248992
#### stand:
# prob
#prob 0.01269835
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 1]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.266"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.146"
checkIC(IC8)
#precision of centering: 6.021289e-10
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 6.021289e-10
Risks(IC8)
#$asCov
# prob
#prob 0.01546324
#
#$asBias
#$asBias$value
#[1] 0.2406692
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.01546324
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.0195648
#
#$asMSE$r
#[1] 0.2661058
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.2661058
plot(IC8)
#-------------------------------------------------------------------------------
## least favorable radius
#-------------------------------------------------------------------------------
## (may take quite some time!)
system.time(r.rho1 <- leastFavorableRadius(L2Fam=B, neighbor=ContNeighborhood(),
risk=asMSE(), rho=0.5))
#current radius: 0.3820278 inefficiency: 1.038184
#current radius: 0.6180722 inefficiency: 1.044504
#current radius: 0.7639556 inefficiency: 1.042432
#current radius: 0.6248193 inefficiency: 1.044545
#current radius: 0.6423138 inefficiency: 1.044575
#current radius: 0.6887768 inefficiency: 1.044172
#current radius: 0.6381993 inefficiency: 1.044581
#current radius: 0.6371031 inefficiency: 1.044575
#current radius: 0.6397095 inefficiency: 1.044575
#current radius: 0.6387762 inefficiency: 1.044584
#current radius: 0.6387355 inefficiency: 1.044584
#current radius: 0.6391327 inefficiency: 1.044573
#current radius: 0.6389123 inefficiency: 1.044572
#current radius: 0.6388282 inefficiency: 1.044571
#current radius: 0.6387762 inefficiency: 1.044584
# user system elapsed
# 125.41 1.20 143.03
r.rho1
#$rho
#[1] 0.5
#
#$leastFavorableRadius
#[1] 0.6387762
#
#$`asMSE-inefficiency`
#[1] 1.044584
system.time(r.rho2 <- leastFavorableRadius(L2Fam=B, neighbor=TotalVarNeighborhood(),
risk=asMSE(), rho=0.5))
#current radius: 0.3820278 inefficiency: 1.040750
#current radius: 0.6180722 inefficiency: 1.028895
#current radius: 0.2361444 inefficiency: 1.043035
#current radius: 0.2225705 inefficiency: 1.042438
#current radius: 0.2881511 inefficiency: 1.043182
#current radius: 0.3240088 inefficiency: 1.043234
#current radius: 0.3524144 inefficiency: 1.042468
#current radius: 0.3103124 inefficiency: 1.043323
#current radius: 0.3081042 inefficiency: 1.043320
#current radius: 0.3106103 inefficiency: 1.043323
#current radius: 0.3105373 inefficiency: 1.043323
#current radius: 0.3104966 inefficiency: 1.043323
#current radius: 0.3105373 inefficiency: 1.043323
# user system elapsed
# 431.94 2.64 483.05
r.rho2
#$rho
#[1] 0.5
#
#$leastFavorableRadius
#[1] 0.3105373
#
#$`asMSE-inefficiency`
#[1] 1.043323
###############################################################################
## k-step (k >= 1) estimation
###############################################################################
## one-step estimation
## 1. generate a contaminated sample
ind <- rbinom(100, size=1, prob=0.05)
x <- rbinom(100, size=25, prob=(1-ind)*0.25 + ind*0.75)
## 2. Kolmogorov(-Smirnov) minimum distance estimator
(est0 <- MDEstimator(x=x, BinomFamily(size=25)))
#Evaluations of Minimum Kolmogorov distance estimate:
#----------------------------------------------------
#An object of class Estimate
#generated by call
# MDEstimator(x = x, ParamFamily = BinomFamily(size = 25))
#samplesize: 100
#estimate:
# prob
#0.2494779
#fixed part of the parameter:
#size
# 25
#Criterion:
#Kolmogorov distance
# 0.05944897
## 3.1. one-step estimation: radius known
## ac) Define infinitesimal robust model
RobB3 <- InfRobModel(center=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(radius=0.5))
## bc) Compute optimally robust IC
IC9 <- optIC(model=RobB3, risk=asMSE())
checkIC(IC9)
#precision of centering: -2.971365e-17
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 2.220446e-16
## cc) Determine 1-step estimate
(est1c <- oneStepEstimator(x, IC=IC9, start=est0))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# oneStepEstimator(x = x, IC = IC9, start = est0)
#samplesize: 100
#estimate:
# prob
# 0.24253514
# (0.00923989)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008537558
#Infos:
# method
#[1,] "oneStepEstimator"
#[2,] "oneStepEstimator"
# message
#[1,] "1-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.06058743
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.249477874158289
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1211749
#### cent: [1] -0.003012083
#### stand:
# prob
#prob 0.01220839
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 1
## instead of ac)-cc) you can also use function roptest
est1c1 <- roptest(x, BinomFamily(size = 25), eps = 0.05, initial.est = est0)
checkIC(pIC(est1c1))
#precision of centering: 4.569788e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
## you can also omit step 2
est1c2 <- roptest(x, BinomFamily(size = 25), eps = 0.05, distance = KolmogorovDist)
checkIC(pIC(est1c2))
#precision of centering: 4.569788e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
## Using Cramer-von-Mises MD estimator (default)
est1c3 <- roptest(x, BinomFamily(size = 25), eps = 0.025)
checkIC(pIC(est1c3))
#precision of centering: -7.770327e-10
#precision of Fisher consistency:
# prob
#prob 2.220446e-16
#maximum deviation
# 7.770327e-10
## comparison of estimates
estimate(est1c)
# prob
#0.2425351
estimate(est1c1)
# prob
#0.2425351
estimate(est1c2)
# prob
#0.2425351
estimate(est1c3)
# prob
#0.2427226
## confidence intervals
confint(est1c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2238655 0.2612048
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC9, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2238655 0.2612048
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, initial.est = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1c2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2238655 0.2612048
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1c3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2251586 0.2602866
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## av) Define infinitesimal robust model
RobB4 <- InfRobModel(center=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(radius=0.25))
## bv) Compute optimally robust IC
IC10 <- optIC(model=RobB4, risk=asMSE())
checkIC(IC10)
## cv) Determine 1-step estimate
(est1v <- oneStepEstimator(x, IC=IC10, start=est0))
## instead of av)-cv) you can also use function roptest
est1v1 <- roptest(x, BinomFamily(size = 25), eps = 0.025, initial.est = est0,
neighbor = TotalVarNeighborhood())
checkIC(pIC(est1v1))
#precision of centering: -9.843168e-15
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
#
9.843168e-15
## you can also omit step 2
est1v2 <- roptest(x, BinomFamily(size = 25), eps = 0.025,
neighbor = TotalVarNeighborhood(), distance = KolmogorovDist)
checkIC(pIC(est1v2))
#precision of centering: -9.843168e-15
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 9.843168e-15
## Using Cramer-von-Mises MD estimator (default)
est1v3 <- roptest(x, BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood())
checkIC(pIC(est1v3))
#precision of centering: 6.67169e-18
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
## comparison of estimates
estimate(est1v)
# prob
#0.2429921
estimate(est1v1)
# prob
#0.2429921
estimate(est1v2)
# prob
#0.2429921
estimate(est1v3)
# prob
#0.2424179
## confidence intervals
confint(est1v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2186101 0.2673741
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC10, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2186101 0.2673741
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, initial.est = est0,
# neighbor = TotalVarNeighborhood())
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1v2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2186101 0.2673741
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood(),
# distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1v3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182528 0.2665830
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood())
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## 3.2. k-step estimation: radius known
IC9 <- optIC(model=RobB3, risk=asMSE())
(est2c <- kStepEstimator(x, IC=IC9, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC9, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.241900493
# (0.009141218)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008356187
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.05992703
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.241953383964985
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1198541
#### cent: [1] -0.004676326
#### stand:
# prob
#prob 0.01194744
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
est2c1 <- roptest(x, BinomFamily(size = 25), eps = 0.05, initial.est = est0, steps = 3L)
checkIC(pIC(est2c1))
#precision of centering: -4.788156e-12
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.788156e-12
est2c2 <- roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3L,
distance = KolmogorovDist)
checkIC(pIC(est2c2))
#precision of centering: -4.788156e-12
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.788156e-12
# Using Cramer-von-Mises MD estimator
est2c3 <- roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3L)
checkIC(pIC(est2c3))
#precision of centering: -6.068637e-12
#precision of Fisher consistency:
# prob
#prob 4.440892e-16
#maximum deviation
# 6.068637e-12
## comparison of estimates
estimate(est2c)
# prob
#0.2419005
estimate(est2c1)
# prob
#0.2419005
estimate(est2c2)
# prob
#0.2419005
estimate(est2c3)
# prob
#0.2418961
## confidence intervals
confint(est2c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234362 0.2603648
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC9, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234362 0.2603648
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, initial.est = est0,
# steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2c2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234362 0.2603648
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, steps = 3L,
# distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2c3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234331 0.2603591
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
IC10 <- optIC(model=RobB4, risk=asMSE())
(est2v <- kStepEstimator(x, IC=IC10, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC10, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24237558
# (0.01195027)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01428089
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.05982794
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.242410409425593
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.25
#
#### clipLo: [1] -0.1148244
#### clipUp: [1] 0.1244874
#### stand:
# prob
#prob 0.01194066
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est2v))
#precision of centering: -2.680039e-18
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
est2v1 <- roptest(x, BinomFamily(size = 25), eps = 0.025, initial.est = est0,
steps = 3L, neighbor = TotalVarNeighborhood())
checkIC(pIC(est2v1))
#precision of centering: 7.299035e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
est2v2 <- roptest(x, BinomFamily(size = 25), eps = 0.025, steps = 3L,
distance = KolmogorovDist, neighbor = TotalVarNeighborhood())
checkIC(pIC(est2v2))
#precision of centering: 7.299035e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
## Using Cramer-von-Mises MD estimator
est2v3 <- roptest(x, BinomFamily(size = 25), eps = 0.025, steps = 3L,
neighbor = TotalVarNeighborhood())
checkIC(pIC(est2v3))
#precision of centering: 3.543879e-18
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
## comparison of estimates
estimate(est2v)
# prob
#0.2423756
estimate(est2v1)
# prob
#0.2423756
estimate(est2v2)
# prob
#0.2423756
estimate(est2v3)
# prob
#0.2423735
## confidence intervals
confint(est2v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182385 0.2665126
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC10, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182385 0.2665126
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, initial.est = est0,
# neighbor = TotalVarNeighborhood(), steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2v2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182385 0.2665126
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood(),
# steps = 3L, distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2v3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182378 0.2665092
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood(),
# steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## 4.1. one-step estimation: radius interval
IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est3c <- oneStepEstimator(x, IC=IC11, start=est0))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# oneStepEstimator(x = x, IC = IC11, start = est0)
#samplesize: 100
#estimate:
# prob
# 0.242379282
# (0.009339235)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.00872213
#Infos:
# method
#[1,] "oneStepEstimator"
#[2,] "oneStepEstimator"
# message
#[1,] "1-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07121993
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.249477874158289
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.6006922
#
#### clip: [1] 0.1185631
#### cent: [1] -0.00453307
#### stand:
# prob
#prob 0.01379441
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, Inf]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.601"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.165"
#steps:
#[1] 1
checkIC(pIC(est3c))
#precision of centering: 9.537084e-18
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 2.220446e-16
IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est3v <- oneStepEstimator(x, IC=IC12, start=est0))
checkIC(pIC(est3v))
## maximum radius for given sample size n: sqrt(n)*0.5
(est3c1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5)
#samplesize: 100
#estimate:
# prob
# 0.241842958
# (0.009234894)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008528328
#Infos:
# method message
#[1,] "roptest" "1-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07044238
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.24311883705201
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5967953
#
#### clip: [1] 0.1180344
#### cent: [1] -0.005169814
#### stand:
# prob
#prob 0.01349046
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 5]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.597"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.159"
#steps:
#[1] 1
checkIC(pIC(est3c1))
#precision of centering: -1.335405e-18
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
(est3v1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5, neighbor = TotalVarNeighborhood()))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood())
#samplesize: 100
#estimate:
# prob
# 0.24231510
# (0.01333293)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01777669
#Infos:
# method message
#[1,] "roptest" "1-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07270451
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.24311883705201
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3104128
#
#### clipLo: [1] -0.1113605
#### clipUp: [1] 0.1228583
#### stand:
# prob
#prob 0.01386519
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 5]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.31"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.166"
#steps:
#[1] 1
checkIC(pIC(est3v1))
#precision of centering: -1.179940e-11
#precision of Fisher consistency:
# prob
#prob 0
#maximum deviation
# 1.179940e-11
## comparison of estimates
estimate(est3c)
# prob
#0.2423793
estimate(est3v)
# prob
#0.2428346
estimate(est3c1)
# prob
#0.2418430
estimate(est3v1)
# prob
#0.2423151
## confidence intervals
confint(est3c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234096 0.261349
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC11, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est3v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2154232 0.2702460
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC12, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est3c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2230924 0.2605935
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est3v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2152137 0.2694165
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood())
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## 4.2. k-step estimation: radius interval
IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est4c <- kStepEstimator(x, IC=IC11, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC11, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24171362
# (0.00922172)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008504013
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07071585
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.241769273363845
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.6006922
#
#### clip: [1] 0.1177239
#### cent: [1] -0.005337597
#### stand:
# prob
#prob 0.01350474
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4c))
#precision of centering: -1.717330e-17
#precision of Fisher consistency:
# prob
#prob 2.220446e-16
#maximum deviation
# 2.220446e-16
IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est4v <- kStepEstimator(x, IC=IC12, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC12, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24224203
# (0.01335018)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01782272
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07293007
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.242288609339931
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3116990
#
#### clipLo: [1] -0.1111696
#### clipUp: [1] 0.1228064
#### stand:
# prob
#prob 0.01387656
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4v))
#precision of centering: -1.454523e-11
#precision of Fisher consistency:
# prob
#prob 4.440892e-16
#maximum deviation
# 1.454523e-11
## maximum radius for given sample size n: sqrt(n)*0.5
(est4c1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.241716814
# (0.009217234)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.00849574
#Infos:
# method message
#[1,] "roptest" "3-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07030448
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.241727413448719
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5967953
#
#### clip: [1] 0.1178033
#### cent: [1] -0.005316906
#### stand:
# prob
#prob 0.01343846
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4c1))
#precision of centering: -1.230939e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
(est4v1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5, neighbor = TotalVarNeighborhood(),
steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood(), steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24224135
# (0.01331745)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01773545
#Infos:
# method message
#[1,] "roptest" "3-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07266429
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.242247135068009
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3104128
#
#### clipLo: [1] -0.1112452
#### clipUp: [1] 0.122844
#### stand:
# prob
#prob 0.01383111
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4v1))
#precision of centering: -1.224149e-11
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.224149e-11
## comparison of estimates
estimate(est4c)
# prob
#0.2417136
estimate(est4v)
# prob
#0.2422420
estimate(est4c1)
# prob
#0.2417168
estimate(est4v1)
# prob
#0.2422414
## confidence intervals
confint(est4c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2229873 0.26044
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC11, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est4v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2151025 0.2693815
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC12, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est4c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2230034 0.2604303
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est4v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2151719 0.2693108
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood(), steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
Try the ROptEst package in your browser
Any scripts or data that you put into this service are public.
ROptEst documentation built on Nov. 17, 2022, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
###############################################################################
## Example: Binomial Family
###############################################################################
require(ROptEst)
options("newDevice"=TRUE)
#-------------------------------------------------------------------------------
## Preparations
#-------------------------------------------------------------------------------
## generates Binomial Family with
## m = 25 and probability of success theta = 0.25
B <- BinomFamily(size = 25, prob = 0.25)
B # show B
#An object of class "BinomFamily"
#### name: Binomial family
#
#### distribution: Distribution Object of Class: Binom
# size: 25
# prob: 0.25
#
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### props:
#[1] "The Binomial family is symmetric with respect to prob = 0.5;"
#[2] "i.e., d(Binom(size, prob))(k)=d(Binom(size,1-prob))(size-k)"
plot(B) # plot of Binom(size = 25, prob = 0.25) and L_2 derivative
checkL2deriv(B)
#precision of centering: -3.103725e-15
#precision of Fisher information:
# prob
#prob -2.842171e-14
#$maximum.deviation
#[1] 2.842171e-14
#-------------------------------------------------------------------------------
## classical optimal IC
#-------------------------------------------------------------------------------
IC0 <- optIC(model = B, risk = asCov())
IC0 # show IC
#An object of class IC
#### name: Classical optimal influence curve for Binomial family
#### L2-differentiable parametric family: Binomial family
#
#### 'Curve': An object of class EuclRandVarList
#Domain: Real Space with dimension 1
#[[1]]
#length of Map: 1
#Range: Real Space with dimension 1
#
#### Infos:
# method message
#[1,] "optIC" "optimal IC in sense of Cramer-Rao bound"
plot(IC0) # plot IC
checkIC(IC0)
#precision of centering: -2.403827e-17
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
Risks(IC0)
#$asCov
# prob
#prob 0.0075
#
#$trAsCov
#[1] 0.0075
#-------------------------------------------------------------------------------
## lower case radius
#-------------------------------------------------------------------------------
lowerCaseRadius(L2Fam = B, neighbor = ContNeighborhood(), risk = asMSE())
#lower case radius
# 1.197809
lowerCaseRadius(L2Fam = B, neighbor = TotalVarNeighborhood(), risk = asMSE())
#lower case radius
# 1.130615
#-------------------------------------------------------------------------------
## L_2 family + infinitesimal neighborhood
#-------------------------------------------------------------------------------
RobB1 <- InfRobModel(center = B, neighbor = ContNeighborhood(radius = 0.5))
RobB1 # show RobB1
#An object of class InfRobModel
####### center: An object of class "BinomFamily"
#### name: Binomial family
#
#### distribution: Distribution Object of Class: Binom
# size: 25
# prob: 0.25
#
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### props:
#[1] "The Binomial family is symmetric with respect to prob = 0.5;"
#[2] "i.e., d(Binom(size, prob))(k)=d(Binom(size,1-prob))(size-k)"
#
####### neighborhood: An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
(RobB2 <- InfRobModel(center = B, neighbor = TotalVarNeighborhood(radius = 0.5)))
#An object of class InfRobModel
####### center: An object of class "BinomFamily"
#### name: Binomial family
#
#### distribution: Distribution Object of Class: Binom
# size: 25
# prob: 0.25
#
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### props:
#[1] "The Binomial family is symmetric with respect to prob = 0.5;"
#[2] "i.e., d(Binom(size, prob))(k)=d(Binom(size,1-prob))(size-k)"
#
####### neighborhood: An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
#-------------------------------------------------------------------------------
## OBRE solution
#-------------------------------------------------------------------------------
system.time(ICA <- optIC(model=RobB1, risk=asAnscombe(),
verbose=TRUE,lower=NULL,upper=10))
#-------------------------------------------------------------------------------
## MSE solution
#-------------------------------------------------------------------------------
system.time(IC1 <- optIC(model=RobB1, risk=asMSE()))
# user system elapsed
# 3.62 0.00 3.62
IC1
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1213661
#### cent: [1] -0.003061948
#### stand:
# prob
#prob 0.01222674
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
checkIC(IC1)
#precision of centering: -3.17283e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
Risks(IC1)
#$asCov
# prob
#prob 0.008544305
#
#$asBias
#$asBias$value
#[1] 0.1213661
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008544305
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01222674
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
getRiskIC(IC1, asBias(), ContNeighborhood()) # standardized bias
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) convex contamination neighborhood"
#
#$asBias$value
#[1] 0.1213661
getRiskIC(IC1, asMSE(), ContNeighborhood(radius = 0.5))
#$asMSE
#$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$neighborhood
#[1] "(uncond.) convex contamination neighborhood with radius 0.5"
#
#$asMSE$radius
#[1] 0.5
(Cov1 <- getRiskIC(IC1, asCov()))
#$asCov
#$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$value
#$asCov$value[[1]]
#[1] 0.008544305
#
#$asCov$value$asCov
#$asCov$value$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$value$asCov$value
# prob
#prob 0.008544305
(mse1 <- getRiskIC(IC1, asMSE(), TotalVarNeighborhood(radius = 0.5)))
#$asMSE
#$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$radius
#[1] 0.5
#
#$asMSE$value
#[1] 0.01222674
(bias1 <- getRiskIC(IC1, asBias(), TotalVarNeighborhood()))
#
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$value
#[1] 0.1213661
## only suboptimal -> ToDo-List
addRisk(IC1) <- list(Cov1, mse1, bias1)
Risks(IC1)
#$asCov
#$asCov[[1]]
#[1] 0.008544305
#
#$asCov$asCov
#$asCov$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value
# prob
#prob 0.008544305
#
#
#$asCov$asCov
#$asCov$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value
#$asCov$asCov$value[[1]]
#[1] 0.008544305
#
#$asCov$asCov$value$asCov
#$asCov$asCov$value$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value$asCov$value
# prob
#prob 0.008544305
#
#
#$asBias
#$asBias$value
#[1] 0.1213661
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#$asBias$asBias
#$asBias$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$asBias$value
#[1] 0.1213661
#
#
#$asBias$asBias
#$asBias$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$asBias$value
#[1] 0.1213661
#
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008544305
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01222674
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
#
#$asMSE$asMSE
#$asMSE$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$asMSE$radius
#[1] 0.5
#
#$asMSE$asMSE$value
#[1] 0.01222674
#
#
#$asMSE$asMSE
#$asMSE$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$asMSE$radius
#[1] 0.5
#
#$asMSE$asMSE$value
#[1] 0.01222674
plot(IC1)
system.time(IC2 <- optIC(model=RobB2, risk=asMSE()))
# user system elapsed
# 9.46 0.02 9.47
IC2
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clipLo: [1] -0.1081481
#### clipUp: [1] 0.1158377
#### stand:
# prob
#prob 0.02208611
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
checkIC(IC2)
#precision of centering: -1.745301e-13
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.745301e-13
Risks(IC2)
#$asCov
# prob
#prob 0.03342380
#
#$asBias
#$asBias$value
#[1] 0.2239858
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.03342380
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.04596621
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
#
getRiskIC(IC2, asMSE(), TotalVarNeighborhood(radius = 0.5))
#$asMSE
#$asMSE$distribution
#[1] "Binom(25, 0.25)"
#
#$asMSE$neighborhood
#[1] "(uncond.) total variation neighborhood with radius 0.5"
#
#$asMSE$radius
#[1] 0.5
#
#$asMSE$value
#[1] 0.04596621
getRiskIC(IC2, asBias(), TotalVarNeighborhood())
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) total variation neighborhood"
#
#$asBias$value
#[1] 0.2239858
getRiskIC(IC2, asBias(), ContNeighborhood())
#$asBias
#$asBias$distribution
#[1] "Binom(25, 0.25)"
#
#$asBias$neighborhood
#[1] "(uncond.) convex contamination neighborhood"
#
#$asBias$value
#[1] 0.2239858
Cov2 <- getRiskIC(IC2, asCov())
addRisk(IC2) <- Cov2
Risks(IC2)
#$asCov
#$asCov[[1]]
#[1] 0.03342380
#
#$asCov$asCov
#$asCov$asCov$distribution
#[1] "Binom(25, 0.25)"
#
#$asCov$asCov$value
# prob
#prob 0.03342380
#
#$asBias
#$asBias$value
#[1] 0.2239858
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.03342380
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.04596621
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
plot(IC2)
#-------------------------------------------------------------------------------
## lower case solutions
#-------------------------------------------------------------------------------
(IC3 <- optIC(model=RobB1, risk=asBias()))
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1098910
#### cent: [1] -1.333333
#### stand:
# [,1]
#[1,] 1
#### lowerCase: [1] -0.3316026
#
#### Infos:
# method message
#[1,] "optIC" "minimum asymptotic bias (lower case) solution"
checkIC(IC3)
#precision of centering: -1.714470e-17
#precision of Fisher consistency:
# prob
#prob 2.220446e-16
#maximum deviation
# 2.220446e-16
Risks(IC3)
#$asBias
#$asBias$value
#[1] 0.1098910
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$asCov
#[1] 0.01011106
#
#$trAsCov
#$trAsCov$value
#[1] 0.01011106
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
#[1] 0.01086581
#
#$asMSE$r
#[1] 0.25
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
plot(IC3)
(IC4 <- optIC(model=RobB2, risk=asBias()))
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
### param: An object of class "ParamFamParameter"
##name: probability of success
#prob: 0.25
##fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clipLo: [1] -0.094766
#### clipUp: [1] 0.1211501
#### stand:
# [,1]
#[1,] 1
#
#### Infos:
# method message
#[1,] "optIC" "minimum asymptotic bias (lower case) solution"
checkIC(IC4)
#precision of centering: -1.076010e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
Risks(IC4)
#$asBias
#$asBias$value
#[1] 0.2159161
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$asCov
#[1] 0.01148091
#
#$trAsCov
#$trAsCov$value
#[1] 0.01148091
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
#[1] 0.01439465
#
#$asMSE$r
#[1] 0.25
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
plot(IC4)
#-------------------------------------------------------------------------------
## Hampel solution
#-------------------------------------------------------------------------------
(IC5 <- optIC(model=RobB1, risk=asHampel(bound=clip(IC1))))
#minimal bound: 0.1098910
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1213661
#### cent: [1] -0.003061961
#### stand:
# prob
#prob 0.01222672
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for 'asHampel' with bound = 0.121"
checkIC(IC5)
#precision of centering: -3.524202e-17
#precision of Fisher consistency:
# prob
#prob -3.56339e-07
#maximum deviation
# 3.56339e-07
Risks(IC5)
#$asCov
# prob
#prob 0.008544296
#
#$asBias
#$asBias$value
#[1] 0.1213661
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008544296
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01222673
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.5
plot(IC5)
(IC6 <- optIC(model=RobB2, risk=asHampel(bound=Risks(IC2)$asBias$value), maxiter = 200))
#minimal bound: 0.2159161
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clipLo: [1] -0.1081481
#### clipUp: [1] 0.1158378
#### stand:
# prob
#prob 0.02208608
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for 'asHampel' with bound = 0.224"
checkIC(IC6)
#precision of centering: -1.745225e-13
#precision of Fisher consistency:
# prob
#prob -1.166564e-07
#maximum deviation
# 1.166564e-07
Risks(IC6)
#$asCov
# prob
#prob 0.03342372
#
#$asBias
#$asBias$value
#[1] 0.2239858
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.03342372
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.04596613
#
#$asMSE$r
#[1] 0.5
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.5
plot(IC6)
#-------------------------------------------------------------------------------
## radius minimax IC
#-------------------------------------------------------------------------------
system.time(IC7 <- radiusMinimaxIC(L2Fam=B, neighbor=ContNeighborhood(),
risk=asMSE(), loRad=0, upRad=1))
# user system elapsed
# 39.39 0.02 39.45
IC7
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3913008
#
#### clip: [1] 0.1269559
#### cent: [1] -0.004305063
#### stand:
# prob
#prob 0.01074017
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 1]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.391"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.103"
checkIC(IC7)
#precision of centering: -3.056847e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
Risks(IC7)
#$asCov
# prob
#prob 0.008272273
#
#$asBias
#$asBias$value
#[1] 0.1269559
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "ContNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.008272273
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.01074017
#
#$asMSE$r
#[1] 0.3913008
#
#$asMSE$at
#An object of class ContNeighborhood
#type: (uncond.) convex contamination neighborhood
#radius: 0.3913008
plot(IC7)
system.time(IC8 <- radiusMinimaxIC(L2Fam=B, neighbor=TotalVarNeighborhood(),
risk=asMSE(), loRad=0, upRad=1))
# user system elapsed
# 163.20 0.12 168.21
IC8
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.25
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.2661058
#
#### clipLo: [1] -0.1157700
#### clipUp: [1] 0.1248992
#### stand:
# prob
#prob 0.01269835
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 1]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.266"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.146"
checkIC(IC8)
#precision of centering: 6.021289e-10
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 6.021289e-10
Risks(IC8)
#$asCov
# prob
#prob 0.01546324
#
#$asBias
#$asBias$value
#[1] 0.2406692
#
#$asBias$biastype
#An object of class symmetricBias
#Slot "name":
#[1] "symmetric Bias"
#
#
#$asBias$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#$asBias$neighbortype
#[1] "TotalVarNeighborhood"
#attr(,"package")
#[1] "RobAStBase"
#
#
#$trAsCov
#$trAsCov$value
# prob
#prob 0.01546324
#
#$trAsCov$normtype
#An object of class NormType
#Slot "name":
#[1] "EuclideanNorm"
#
#Slot "fct":
#function (x)
#{
# if (is.vector(x))
# return(abs(x))
# else return(sqrt(colSums(x^2)))
#}
#<environment: namespace:distrMod>
#
#
#
#$asMSE
#$asMSE$value
# prob
#prob 0.0195648
#
#$asMSE$r
#[1] 0.2661058
#
#$asMSE$at
#An object of class TotalVarNeighborhood
#type: (uncond.) total variation neighborhood
#radius: 0.2661058
plot(IC8)
#-------------------------------------------------------------------------------
## least favorable radius
#-------------------------------------------------------------------------------
## (may take quite some time!)
system.time(r.rho1 <- leastFavorableRadius(L2Fam=B, neighbor=ContNeighborhood(),
risk=asMSE(), rho=0.5))
#current radius: 0.3820278 inefficiency: 1.038184
#current radius: 0.6180722 inefficiency: 1.044504
#current radius: 0.7639556 inefficiency: 1.042432
#current radius: 0.6248193 inefficiency: 1.044545
#current radius: 0.6423138 inefficiency: 1.044575
#current radius: 0.6887768 inefficiency: 1.044172
#current radius: 0.6381993 inefficiency: 1.044581
#current radius: 0.6371031 inefficiency: 1.044575
#current radius: 0.6397095 inefficiency: 1.044575
#current radius: 0.6387762 inefficiency: 1.044584
#current radius: 0.6387355 inefficiency: 1.044584
#current radius: 0.6391327 inefficiency: 1.044573
#current radius: 0.6389123 inefficiency: 1.044572
#current radius: 0.6388282 inefficiency: 1.044571
#current radius: 0.6387762 inefficiency: 1.044584
# user system elapsed
# 125.41 1.20 143.03
r.rho1
#$rho
#[1] 0.5
#
#$leastFavorableRadius
#[1] 0.6387762
#
#$`asMSE-inefficiency`
#[1] 1.044584
system.time(r.rho2 <- leastFavorableRadius(L2Fam=B, neighbor=TotalVarNeighborhood(),
risk=asMSE(), rho=0.5))
#current radius: 0.3820278 inefficiency: 1.040750
#current radius: 0.6180722 inefficiency: 1.028895
#current radius: 0.2361444 inefficiency: 1.043035
#current radius: 0.2225705 inefficiency: 1.042438
#current radius: 0.2881511 inefficiency: 1.043182
#current radius: 0.3240088 inefficiency: 1.043234
#current radius: 0.3524144 inefficiency: 1.042468
#current radius: 0.3103124 inefficiency: 1.043323
#current radius: 0.3081042 inefficiency: 1.043320
#current radius: 0.3106103 inefficiency: 1.043323
#current radius: 0.3105373 inefficiency: 1.043323
#current radius: 0.3104966 inefficiency: 1.043323
#current radius: 0.3105373 inefficiency: 1.043323
# user system elapsed
# 431.94 2.64 483.05
r.rho2
#$rho
#[1] 0.5
#
#$leastFavorableRadius
#[1] 0.3105373
#
#$`asMSE-inefficiency`
#[1] 1.043323
###############################################################################
## k-step (k >= 1) estimation
###############################################################################
## one-step estimation
## 1. generate a contaminated sample
ind <- rbinom(100, size=1, prob=0.05)
x <- rbinom(100, size=25, prob=(1-ind)*0.25 + ind*0.75)
## 2. Kolmogorov(-Smirnov) minimum distance estimator
(est0 <- MDEstimator(x=x, BinomFamily(size=25)))
#Evaluations of Minimum Kolmogorov distance estimate:
#----------------------------------------------------
#An object of class Estimate
#generated by call
# MDEstimator(x = x, ParamFamily = BinomFamily(size = 25))
#samplesize: 100
#estimate:
# prob
#0.2494779
#fixed part of the parameter:
#size
# 25
#Criterion:
#Kolmogorov distance
# 0.05944897
## 3.1. one-step estimation: radius known
## ac) Define infinitesimal robust model
RobB3 <- InfRobModel(center=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(radius=0.5))
## bc) Compute optimally robust IC
IC9 <- optIC(model=RobB3, risk=asMSE())
checkIC(IC9)
#precision of centering: -2.971365e-17
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 2.220446e-16
## cc) Determine 1-step estimate
(est1c <- oneStepEstimator(x, IC=IC9, start=est0))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# oneStepEstimator(x = x, IC = IC9, start = est0)
#samplesize: 100
#estimate:
# prob
# 0.24253514
# (0.00923989)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008537558
#Infos:
# method
#[1,] "oneStepEstimator"
#[2,] "oneStepEstimator"
# message
#[1,] "1-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.06058743
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.249477874158289
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1211749
#### cent: [1] -0.003012083
#### stand:
# prob
#prob 0.01220839
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 1
## instead of ac)-cc) you can also use function roptest
est1c1 <- roptest(x, BinomFamily(size = 25), eps = 0.05, initial.est = est0)
checkIC(pIC(est1c1))
#precision of centering: 4.569788e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
## you can also omit step 2
est1c2 <- roptest(x, BinomFamily(size = 25), eps = 0.05, distance = KolmogorovDist)
checkIC(pIC(est1c2))
#precision of centering: 4.569788e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
## Using Cramer-von-Mises MD estimator (default)
est1c3 <- roptest(x, BinomFamily(size = 25), eps = 0.025)
checkIC(pIC(est1c3))
#precision of centering: -7.770327e-10
#precision of Fisher consistency:
# prob
#prob 2.220446e-16
#maximum deviation
# 7.770327e-10
## comparison of estimates
estimate(est1c)
# prob
#0.2425351
estimate(est1c1)
# prob
#0.2425351
estimate(est1c2)
# prob
#0.2425351
estimate(est1c3)
# prob
#0.2427226
## confidence intervals
confint(est1c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2238655 0.2612048
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC9, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2238655 0.2612048
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, initial.est = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1c2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2238655 0.2612048
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1c3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2251586 0.2602866
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## av) Define infinitesimal robust model
RobB4 <- InfRobModel(center=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(radius=0.25))
## bv) Compute optimally robust IC
IC10 <- optIC(model=RobB4, risk=asMSE())
checkIC(IC10)
## cv) Determine 1-step estimate
(est1v <- oneStepEstimator(x, IC=IC10, start=est0))
## instead of av)-cv) you can also use function roptest
est1v1 <- roptest(x, BinomFamily(size = 25), eps = 0.025, initial.est = est0,
neighbor = TotalVarNeighborhood())
checkIC(pIC(est1v1))
#precision of centering: -9.843168e-15
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
#
9.843168e-15
## you can also omit step 2
est1v2 <- roptest(x, BinomFamily(size = 25), eps = 0.025,
neighbor = TotalVarNeighborhood(), distance = KolmogorovDist)
checkIC(pIC(est1v2))
#precision of centering: -9.843168e-15
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 9.843168e-15
## Using Cramer-von-Mises MD estimator (default)
est1v3 <- roptest(x, BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood())
checkIC(pIC(est1v3))
#precision of centering: 6.67169e-18
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
## comparison of estimates
estimate(est1v)
# prob
#0.2429921
estimate(est1v1)
# prob
#0.2429921
estimate(est1v2)
# prob
#0.2429921
estimate(est1v3)
# prob
#0.2424179
## confidence intervals
confint(est1v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2186101 0.2673741
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC10, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2186101 0.2673741
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, initial.est = est0,
# neighbor = TotalVarNeighborhood())
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1v2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2186101 0.2673741
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood(),
# distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est1v3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182528 0.2665830
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood())
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## 3.2. k-step estimation: radius known
IC9 <- optIC(model=RobB3, risk=asMSE())
(est2c <- kStepEstimator(x, IC=IC9, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC9, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.241900493
# (0.009141218)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008356187
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.05992703
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.241953383964985
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5
#
#### clip: [1] 0.1198541
#### cent: [1] -0.004676326
#### stand:
# prob
#prob 0.01194744
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
est2c1 <- roptest(x, BinomFamily(size = 25), eps = 0.05, initial.est = est0, steps = 3L)
checkIC(pIC(est2c1))
#precision of centering: -4.788156e-12
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.788156e-12
est2c2 <- roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3L,
distance = KolmogorovDist)
checkIC(pIC(est2c2))
#precision of centering: -4.788156e-12
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.788156e-12
# Using Cramer-von-Mises MD estimator
est2c3 <- roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3L)
checkIC(pIC(est2c3))
#precision of centering: -6.068637e-12
#precision of Fisher consistency:
# prob
#prob 4.440892e-16
#maximum deviation
# 6.068637e-12
## comparison of estimates
estimate(est2c)
# prob
#0.2419005
estimate(est2c1)
# prob
#0.2419005
estimate(est2c2)
# prob
#0.2419005
estimate(est2c3)
# prob
#0.2418961
## confidence intervals
confint(est2c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234362 0.2603648
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC9, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234362 0.2603648
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, initial.est = est0,
# steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2c2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234362 0.2603648
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, steps = 3L,
# distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2c3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234331 0.2603591
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.05, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
IC10 <- optIC(model=RobB4, risk=asMSE())
(est2v <- kStepEstimator(x, IC=IC10, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC10, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24237558
# (0.01195027)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01428089
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.05982794
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.242410409425593
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.25
#
#### clipLo: [1] -0.1148244
#### clipUp: [1] 0.1244874
#### stand:
# prob
#prob 0.01194066
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est2v))
#precision of centering: -2.680039e-18
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
est2v1 <- roptest(x, BinomFamily(size = 25), eps = 0.025, initial.est = est0,
steps = 3L, neighbor = TotalVarNeighborhood())
checkIC(pIC(est2v1))
#precision of centering: 7.299035e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
est2v2 <- roptest(x, BinomFamily(size = 25), eps = 0.025, steps = 3L,
distance = KolmogorovDist, neighbor = TotalVarNeighborhood())
checkIC(pIC(est2v2))
#precision of centering: 7.299035e-18
#precision of Fisher consistency:
# prob
#prob -4.440892e-16
#maximum deviation
# 4.440892e-16
## Using Cramer-von-Mises MD estimator
est2v3 <- roptest(x, BinomFamily(size = 25), eps = 0.025, steps = 3L,
neighbor = TotalVarNeighborhood())
checkIC(pIC(est2v3))
#precision of centering: 3.543879e-18
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
## comparison of estimates
estimate(est2v)
# prob
#0.2423756
estimate(est2v1)
# prob
#0.2423756
estimate(est2v2)
# prob
#0.2423756
estimate(est2v3)
# prob
#0.2423735
## confidence intervals
confint(est2v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182385 0.2665126
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC10, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182385 0.2665126
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, initial.est = est0,
# neighbor = TotalVarNeighborhood(), steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2v2, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182385 0.2665126
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood(),
# steps = 3L, distance = KolmogorovDist)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est2v3, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2182378 0.2665092
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps = 0.025, neighbor = TotalVarNeighborhood(),
# steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## 4.1. one-step estimation: radius interval
IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est3c <- oneStepEstimator(x, IC=IC11, start=est0))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# oneStepEstimator(x = x, IC = IC11, start = est0)
#samplesize: 100
#estimate:
# prob
# 0.242379282
# (0.009339235)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.00872213
#Infos:
# method
#[1,] "oneStepEstimator"
#[2,] "oneStepEstimator"
# message
#[1,] "1-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07121993
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.249477874158289
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.6006922
#
#### clip: [1] 0.1185631
#### cent: [1] -0.00453307
#### stand:
# prob
#prob 0.01379441
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, Inf]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.601"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.165"
#steps:
#[1] 1
checkIC(pIC(est3c))
#precision of centering: 9.537084e-18
#precision of Fisher consistency:
# prob
#prob -2.220446e-16
#maximum deviation
# 2.220446e-16
IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est3v <- oneStepEstimator(x, IC=IC12, start=est0))
checkIC(pIC(est3v))
## maximum radius for given sample size n: sqrt(n)*0.5
(est3c1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5)
#samplesize: 100
#estimate:
# prob
# 0.241842958
# (0.009234894)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008528328
#Infos:
# method message
#[1,] "roptest" "1-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07044238
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.24311883705201
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5967953
#
#### clip: [1] 0.1180344
#### cent: [1] -0.005169814
#### stand:
# prob
#prob 0.01349046
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 5]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.597"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.159"
#steps:
#[1] 1
checkIC(pIC(est3c1))
#precision of centering: -1.335405e-18
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.110223e-16
(est3v1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5, neighbor = TotalVarNeighborhood()))
#Evaluations of 1-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood())
#samplesize: 100
#estimate:
# prob
# 0.24231510
# (0.01333293)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01777669
#Infos:
# method message
#[1,] "roptest" "1-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07270451
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.24311883705201
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3104128
#
#### clipLo: [1] -0.1113605
#### clipUp: [1] 0.1228583
#### stand:
# prob
#prob 0.01386519
#
#### Infos:
# method message
#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0, 5]"
#[2,] "radiusMinimaxIC" "least favorable radius: 0.31"
#[3,] "radiusMinimaxIC" "maximum asMSE-inefficiency: 1.166"
#steps:
#[1] 1
checkIC(pIC(est3v1))
#precision of centering: -1.179940e-11
#precision of Fisher consistency:
# prob
#prob 0
#maximum deviation
# 1.179940e-11
## comparison of estimates
estimate(est3c)
# prob
#0.2423793
estimate(est3v)
# prob
#0.2428346
estimate(est3c1)
# prob
#0.2418430
estimate(est3v1)
# prob
#0.2423151
## confidence intervals
confint(est3c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2234096 0.261349
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC11, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est3v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2154232 0.2702460
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#oneStepEstimator(x = x, IC = IC12, start = est0)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est3c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2230924 0.2605935
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est3v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2152137 0.2694165
#Type of estimator: 1-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood())
#Fixed part of the parameter at which estimate was produced:
#size
# 25
## 4.2. k-step estimation: radius interval
IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est4c <- kStepEstimator(x, IC=IC11, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC11, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24171362
# (0.00922172)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.008504013
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07071585
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.241769273363845
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.6006922
#
#### clip: [1] 0.1177239
#### cent: [1] -0.005337597
#### stand:
# prob
#prob 0.01350474
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4c))
#precision of centering: -1.717330e-17
#precision of Fisher consistency:
# prob
#prob 2.220446e-16
#maximum deviation
# 2.220446e-16
IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
(est4v <- kStepEstimator(x, IC=IC12, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# kStepEstimator(x = x, IC = IC12, start = est0, steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24224203
# (0.01335018)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01782272
#Infos:
# method
#[1,] "kStepEstimator"
#[2,] "kStepEstimator"
# message
#[1,] "3-step estimate for Binomial family"
#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07293007
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.242288609339931
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3116990
#
#### clipLo: [1] -0.1111696
#### clipUp: [1] 0.1228064
#### stand:
# prob
#prob 0.01387656
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4v))
#precision of centering: -1.454523e-11
#precision of Fisher consistency:
# prob
#prob 4.440892e-16
#maximum deviation
# 1.454523e-11
## maximum radius for given sample size n: sqrt(n)*0.5
(est4c1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.241716814
# (0.009217234)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.00849574
#Infos:
# method message
#[1,] "roptest" "3-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07030448
#(partial) influence curve:
#An object of class ContIC
#### name: IC of contamination type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.241727413448719
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.5967953
#
#### clip: [1] 0.1178033
#### cent: [1] -0.005316906
#### stand:
# prob
#prob 0.01343846
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4c1))
#precision of centering: -1.230939e-17
#precision of Fisher consistency:
# prob
#prob -3.330669e-16
#maximum deviation
# 3.330669e-16
(est4v1 <- roptest(x, BinomFamily(size = 25), eps.upper = 0.5, neighbor = TotalVarNeighborhood(),
steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
#An object of class Estimate
#generated by call
# roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood(), steps = 3L)
#samplesize: 100
#estimate:
# prob
# 0.24224135
# (0.01331745)
#fixed part of the parameter:
#size
# 25
#asymptotic (co)variance (multiplied with samplesize):
#[1] 0.01773545
#Infos:
# method message
#[1,] "roptest" "3-step estimate for Binomial family"
#[2,] "roptest" "computation of IC, asvar and asbias via useLast = FALSE"
#asymptotic bias:
#[1] 0.07266429
#(partial) influence curve:
#An object of class TotalVarIC
#### name: IC of total variation type
#
#### L2-differentiable parametric family: Binomial family
#### param: An object of class "ParamFamParameter"
#name: probability of success
#prob: 0.242247135068009
#fixed part of param.:
# size: 25
#trafo:
# prob
#prob 1
#
#### neighborhood radius: 0.3104128
#
#### clipLo: [1] -0.1112452
#### clipUp: [1] 0.122844
#### stand:
# prob
#prob 0.01383111
#
#### Infos:
# method message
#[1,] "optIC" "optimally robust IC for asMSE"
#steps:
#[1] 3
checkIC(pIC(est4v1))
#precision of centering: -1.224149e-11
#precision of Fisher consistency:
# prob
#prob -1.110223e-16
#maximum deviation
# 1.224149e-11
## comparison of estimates
estimate(est4c)
# prob
#0.2417136
estimate(est4v)
# prob
#0.2422420
estimate(est4c1)
# prob
#0.2417168
estimate(est4v1)
# prob
#0.2422414
## confidence intervals
confint(est4c, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2229873 0.26044
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC11, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est4v, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2151025 0.2693815
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#kStepEstimator(x = x, IC = IC12, start = est0, steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est4c1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2230034 0.2604303
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
confint(est4v1, symmetricBias())
#A[n] asymptotic (LAN-based), uniform (bias-aware)
# confidence interval:
#for symmetric Bias
# 2.5 % 97.5 %
#prob 0.2151719 0.2693108
#Type of estimator: 3-step estimate
#samplesize: 100
#Call by which estimate was produced:
#roptest(x = x, L2Fam = BinomFamily(size = 25), eps.upper = 0.5,
# neighbor = TotalVarNeighborhood(), steps = 3L)
#Fixed part of the parameter at which estimate was produced:
#size
# 25
Try the ROptEst package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.