[Robast-commits] r362 - in branches/robast-0.7/pkg: ROptEst/chm ROptEst/inst/scripts RobAStBase/R RobAStBase/chm RobAStBase/man RobLox/R RobLox/chm
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 1 06:35:14 CEST 2009
Author: ruckdeschel
Date: 2009-09-01 06:35:13 +0200 (Tue, 01 Sep 2009)
New Revision: 362
Modified:
branches/robast-0.7/pkg/ROptEst/chm/ROptEst.chm
branches/robast-0.7/pkg/ROptEst/inst/scripts/BinomialModel.R
branches/robast-0.7/pkg/ROptEst/inst/scripts/PoissonModel.R
branches/robast-0.7/pkg/RobAStBase/R/AllClass.R
branches/robast-0.7/pkg/RobAStBase/R/kStepEstimator.R
branches/robast-0.7/pkg/RobAStBase/R/oneStepEstimator.R
branches/robast-0.7/pkg/RobAStBase/chm/RobAStBase.chm
branches/robast-0.7/pkg/RobAStBase/chm/kStepEstimator.html
branches/robast-0.7/pkg/RobAStBase/chm/oneStepEstimator.html
branches/robast-0.7/pkg/RobAStBase/man/kStepEstimator.Rd
branches/robast-0.7/pkg/RobAStBase/man/oneStepEstimator.Rd
branches/robast-0.7/pkg/RobLox/R/roblox.R
branches/robast-0.7/pkg/RobLox/R/rowRoblox.R
branches/robast-0.7/pkg/RobLox/chm/RobLox.chm
Log:
RobLox: gains new slots of kStepEstimator
ROptEst: scripts in part now with R output in comment
RobAStBase: + StartClass now also contains "matrix"
+ argument start in kStepEstimator and oneStepEstimator is NULL by default (and then replaced by L2Fam at startPar)
Modified: branches/robast-0.7/pkg/ROptEst/chm/ROptEst.chm
===================================================================
(Binary files differ)
Modified: branches/robast-0.7/pkg/ROptEst/inst/scripts/BinomialModel.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/inst/scripts/BinomialModel.R 2009-08-31 19:53:18 UTC (rev 361)
+++ branches/robast-0.7/pkg/ROptEst/inst/scripts/BinomialModel.R 2009-09-01 04:35:13 UTC (rev 362)
@@ -1,122 +1,1490 @@
###############################################################################
## Example: Binomial Family
###############################################################################
+
require(ROptEst)
options("newDevice"=TRUE)
-## generates Binomial Family with
+#-------------------------------------------------------------------------------
+## 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
+
+#-------------------------------------------------------------------------------
## 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
+
+###ERROR
(mse1 <- getRiskIC(IC1, asMSE(), TotalVarNeighborhood(radius = 0.5)))
-(bias1 <- getRiskIC(IC1, asBias(), TotalVarNeighborhood()))
+
+#(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(),
+#-------------------------------------------------------------------------------
+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(),
+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
+# 437.47 0.09 438.58
+
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
+#2211.92 1.13 2235.70
+
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)
+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)),
@@ -124,31 +1492,169 @@
## 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)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 362
More information about the Robast-commits
mailing list