[Robast-commits] r407 - branches/robast-0.8/pkg/ROptEst/inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 28 22:35:02 CEST 2010
Author: ruckdeschel
Date: 2010-05-28 22:35:02 +0200 (Fri, 28 May 2010)
New Revision: 407
Added:
branches/robast-0.8/pkg/ROptEst/inst/scripts/NbinomModel.R
Log:
uploaded Script for NbinomFamily
Added: branches/robast-0.8/pkg/ROptEst/inst/scripts/NbinomModel.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/inst/scripts/NbinomModel.R (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/inst/scripts/NbinomModel.R 2010-05-28 20:35:02 UTC (rev 407)
@@ -0,0 +1,2812 @@
+###############################################################################
+## Example: Neg Binomial Family
+###############################################################################
+
+require(ROptEst)
+options("newDevice"=TRUE)
+
+#-------------------------------------------------------------------------------
+## Preparations
+#-------------------------------------------------------------------------------
+## generates Neg.Binomial Family with
+## m = 25 and probability of success theta = 0.25
+N <- NbinomFamily(size = 25, prob = 0.25)
+N # show N
+
+#An object of class "NbinomFamily"
+#### name: Negative Binomial family
+#
+#### distribution: Distribution Object of Class: Nbinom
+# 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] ""
+#
+#plot(N) # plot of Nbinom(size = 25, prob = 0.25) and L_2 derivative
+#checkL2deriv(N)
+#
+##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 = N, risk = asCov())
+IC0 # show IC
+
+#An object of class “IC”
+#### name: Classical optimal influence curve for Negative Binomial family
+#### L2-differentiable parametric family: Negative 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: 7.796853e-15
+#precision of Fisher consistency:
+# prob
+#prob -2.166600e-12
+#maximum deviation
+# 2.166600e-12
+
+Risks(IC0)
+
+#$asCov
+# prob
+#prob 0.001875
+#
+#$trAsCov
+#[1] 0.001875
+
+#-------------------------------------------------------------------------------
+## lower case radius
+#-------------------------------------------------------------------------------
+lowerCaseRadius(L2Fam = N, neighbor = ContNeighborhood(), risk = asMSE())
+
+#lower case radius
+# 4.153322
+
+lowerCaseRadius(L2Fam = N, neighbor = TotalVarNeighborhood(), risk = asMSE())
+
+#lower case radius
+# 1.840705
+
+#-------------------------------------------------------------------------------
+## L_2 family + infinitesimal neighborhood
+#-------------------------------------------------------------------------------
+RobN1 <- InfRobModel(center = N, neighbor = ContNeighborhood(radius = 0.5))
+RobN1 # show RobN1
+
+#An object of class “InfRobModel”
+####### center: An object of class "NbinomFamily"
+#### name: Negative Binomial family
+#
+#### distribution: Distribution Object of Class: Nbinom
+# 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] ""
+#
+####### neighborhood: An object of class “ContNeighborhood”
+#type: (uncond.) convex contamination neighborhood
+#radius: 0.5
+
+(RobN2 <- InfRobModel(center = N, neighbor = TotalVarNeighborhood(radius = 0.5)))
+
+#An object of class “InfRobModel”
+####### center: An object of class "NbinomFamily"
+#### name: Negative Binomial family
+#
+#### distribution: Distribution Object of Class: Nbinom
+# 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] ""
+#
+####### neighborhood: An object of class “TotalVarNeighborhood”
+#type: (uncond.) total variation neighborhood
+#radius: 0.5
+
+#-------------------------------------------------------------------------------
+## MSE solution
+#-------------------------------------------------------------------------------
+system.time(IC1 <- optIC(model=RobN1, risk=asMSE()))
+
+# user system elapsed
+# 10.53 0.02 11.21
+
+IC1
+
+#An object of class “ContIC”
+#### name: IC of contamination type
+#
+#### L2-differentiable parametric family: Negative 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.06143505
+#### cent: [1] 0.003717104
+#### stand:
+# prob
+#prob 0.00310076
+#
+#### Infos:
+# method message
+#[1,] "optIC" "optimally robust IC for ‘asMSE’"
+
+checkIC(IC1)
+
+#precision of centering: 4.297326e-14
+#precision of Fisher consistency:
+# prob
+#prob 2.555733e-13
+#maximum deviation
+# 2.555733e-13
+
+Risks(IC1)
+
+#$asCov
+# prob
+#prob 0.002157193
+#
+#$asBias
+#$asBias$value
+#[1] 0.06143505
+#
+#$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.002157193
+#
+#$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.00310076
+#
+#$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] "Nbinom(25, 0.25)"
+#
+#$asBias$neighborhood
+#[1] "(uncond.) convex contamination neighborhood"
+#
+#$asBias$value
+#[1] 0.06143505
+
+
+getRiskIC(IC1, asMSE(), ContNeighborhood(radius = 0.5))
+
+#$asMSE
+#$asMSE$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asMSE$neighborhood
+#[1] "(uncond.) convex contamination neighborhood with radius 0.5"
+#
+#$asMSE$radius
+#[1] 0.5
+#
+#$asMSE$value
+#[1] 0.00310076
+
+
+(Cov1 <- getRiskIC(IC1, asCov()))
+
+#$asCov
+#$asCov$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asCov$value
+# prob
+#prob 0.002157193
+
+
+
+(mse1 <- getRiskIC(IC1, asMSE(), TotalVarNeighborhood(radius = 0.5)))
+
+#$asMSE
+#$asMSE$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asMSE$neighborhood
+#[1] "(uncond.) total variation neighborhood with radius 0.5"
+#
+#$asMSE$radius
+#[1] 0.5
+#
+#$asMSE$value
+#[1] 0.00310076
+
+(bias1 <- getRiskIC(IC1, asBias(), TotalVarNeighborhood()))
+
+#$asBias
+#$asBias$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asBias$neighborhood
+#[1] "(uncond.) total variation neighborhood"
+#
+#$asBias$value
+#[1] 0.06143505
+
+
+## only suboptimal -> ToDo-List
+addRisk(IC1) <- list(Cov1, mse1, bias1)
+Risks(IC1)
+
+#$asCov
+#$asCov[[1]]
+#[1] 0.002157193
+#
+#$asCov$asCov
+#$asCov$asCov$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asCov$asCov$value
+# prob
+#prob 0.002157193
+#
+#
+#
+#$asBias
+#$asBias$value
+#[1] 0.06143505
+#
+#$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] "Nbinom(25, 0.25)"
+#
+#$asBias$asBias$neighborhood
+#[1] "(uncond.) total variation neighborhood"
+#
+#$asBias$asBias$value
+#[1] 0.06143505
+#
+#
+#
+#$trAsCov
+#$trAsCov$value
+# prob
+#prob 0.002157193
+#
+#$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.00310076
+#
+#$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] "Nbinom(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.00310076
+
+
+plot(IC1)
+
+system.time(IC2 <- optIC(model=RobN2, risk=asMSE()))
+
+# user system elapsed
+# 75.57 0.22 81.59
+
+IC2
+
+#An object of class “TotalVarIC”
+#### name: IC of total variation type
+#
+#### L2-differentiable parametric family: Negative 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.06032159
+#### clipUp: [1] 0.052032
+#### stand:
+# prob
+#prob 0.005585238
+#
+#### Infos:
+# method message
+#[1,] "optIC" "optimally robust IC for ‘asMSE’"
+
+checkIC(IC2)
+
+#precision of centering: 9.638494e-16
+#precision of Fisher consistency:
+# prob
+#prob 2.164935e-13
+#maximum deviation
+# 2.164935e-13
+
+Risks(IC2)
+
+#$asCov
+# prob
+#prob 0.008255488
+#
+#$asBias
+#$asBias$value
+#[1] 0.1123536
+#
+#$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.008255488
+#
+#$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.01141132
+#
+#$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] "Nbinom(25, 0.25)"
+#
+#$asMSE$neighborhood
+#[1] "(uncond.) total variation neighborhood with radius 0.5"
+#
+#$asMSE$radius
+#[1] 0.5
+#
+#$asMSE$value
+#[1] 0.01141132
+
+getRiskIC(IC2, asBias(), TotalVarNeighborhood())
+
+#$asBias
+#$asBias$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asBias$neighborhood
+#[1] "(uncond.) total variation neighborhood"
+#
+#$asBias$value
+#[1] 0.1123536
+
+getRiskIC(IC2, asBias(), ContNeighborhood())
+
+#$asBias
+#$asBias$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asBias$neighborhood
+#[1] "(uncond.) convex contamination neighborhood"
+#
+#$asBias$value
+#[1] 0.1123536
+
+Cov2 <- getRiskIC(IC2, asCov())
+addRisk(IC2) <- Cov2
+Risks(IC2)
+
+#$asCov
+#$asCov[[1]]
+#[1] 0.008255488
+#
+#$asCov$asCov
+#$asCov$asCov$distribution
+#[1] "Nbinom(25, 0.25)"
+#
+#$asCov$asCov$value
+# prob
+#prob 0.008255488
+#
+#
+#
+#$asBias
+#$asBias$value
+#[1] 0.1123536
+#
+#$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.008255488
+#
+#$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.01141132
+#
+#$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=RobN1, risk=asBias()))
+
+
+#An object of class “ContIC”
+#### name: IC of contamination type
+#
+#### L2-differentiable parametric family: Negative 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.05458835
+#### cent: [1] 1.333333
+#### stand:
+# [,1]
+#[1,] 1
+#### lowerCase: [1] -0.3268154
+#
+#### Infos:
+# method message
+#[1,] "optIC" "minimum asymptotic bias (lower case) solution"
+
+
+checkIC(IC3)
+
+#precision of centering: 7.555495e-16
+#precision of Fisher consistency:
+# prob
+#prob 8.881784e-16
+#maximum deviation
+# 8.881784e-16
+
+Risks(IC3)
+
+#$asBias
+#$asBias$value
+#[1] 0.05458835
+#
+#$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.002918187
+#
+#$trAsCov
+#$trAsCov$value
+#[1] 0.002918187
+#
+#$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.00310443
+#
+#$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=RobN2, risk=asBias()))
+
+#An object of class “TotalVarIC”
+#### name: IC of total variation type
+#
+#### L2-differentiable parametric family: Negative 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.0574604
+#### clipUp: [1] 0.05147243
+#### stand:
+# [,1]
+#[1,] 1
+#
+#### Infos:
+# method message
+#[1,] "optIC" "minimum asymptotic bias (lower case) solution"
+
+checkIC(IC4)
+
+#precision of centering: 1.111799e-15
+#precision of Fisher consistency:
+# prob
+#prob 2.14051e-13
+#maximum deviation
+# 2.14051e-13
+
+Risks(IC4)
+
+#$asBias
+#$asBias$value
+#[1] 0.1089328
+#
+#$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.002889749
+#
+#$trAsCov
+#$trAsCov$value
+#[1] 0.002889749
+#
+#$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.003631397
+#
+#$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=RobN1, risk=asHampel(bound=clip(IC1))))
+
+#minimal bound: 0.05458835
+#An object of class “ContIC”
+#### name: IC of contamination type
+#
+#### L2-differentiable parametric family: Negative 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.06143505
+#### cent: [1] 0.003717089
+#### stand:
+# prob
+#prob 0.003100752
+#
+#### Infos:
+# method message
+#[1,] "optIC" "optimally robust IC for 'asHampel' with bound = 0.061"
+
+
+checkIC(IC5)
+
+#precision of centering: 4.289037e-14
+#precision of Fisher consistency:
+# prob
+#prob -5.088488e-07
+#maximum deviation
+# 5.088488e-07
+
+Risks(IC5)
+
+#$asCov
+# prob
+#prob 0.002157190
+#
+#$asBias
+#$asBias$value
+#[1] 0.06143505
+#
+#$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.002157190
+#
+#$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.003100757
+#
+#$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=RobN2, risk=asHampel(bound=Risks(IC2)$asBias$value), maxiter = 200))
+
+#minimal bound: 0.1089328
+#maximum iterations reached!
+# achieved precision: 5.583403e-07
+#An object of class “TotalVarIC”
+#### name: IC of total variation type
+#
+#### L2-differentiable parametric family: Negative 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.06032159
+#### clipUp: [1] 0.052032
+#### stand:
+# prob
+#prob 0.005585234
+#
+#### Infos:
+# method message
+#[1,] "optIC" "optimally robust IC for 'asHampel' with bound = 0.112"
+
+
+checkIC(IC6)
+
+#precision of centering: 9.55774e-16
+#precision of Fisher consistency:
+# prob
+#prob -4.638466e-08
+#maximum deviation
+# 4.638466e-08
+
+
+Risks(IC6)
+
+#$asCov
+# prob
+#prob 0.008255479
+#
+#$asBias
+#$asBias$value
+#[1] 0.1123536
+#
+#$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.008255479
+#
+#$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.01141131
+#
+#$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=N, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0.01, upRad=3.9))
+# user system elapsed
+# 33.26 0.02 33.64
+
+IC7
+
+#An object of class “ContIC”
+#### name: IC of contamination type
+#
+#### L2-differentiable parametric family: Negative 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.5814856
+#
+#### clip: [1] 0.05991006
+#### cent: [1] 0.004363098
+#### stand:
+# prob
+#prob 0.003424638
+#
+#### Infos:
+# method message
+#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0.01, 3.9]"
+#[2,] "radiusMinimaxIC" "least favorable radius: 0.581"
+#[3,] "radiusMinimaxIC" "maximum ‘asMSE’-inefficiency: 1.177"
+
+checkIC(IC7)
+
+#precision of centering: 5.901277e-12
+#precision of Fisher consistency:
+# prob
+#prob -1.295138e-09
+#maximum deviation
+# 1.295138e-09
+
+Risks(IC7)
+
+#$asCov
+# prob
+#prob 0.002211033
+#
+#$asBias
+#$asBias$value
+#[1] 0.05991006
+#
+#$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.002211033
+#
+#$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.003424638
+#
+#$asMSE$r
+#[1] 0.5814856
+#
+#$asMSE$at
+#An object of class “ContNeighborhood”
+#type: (uncond.) convex contamination neighborhood
+#radius: 0.5814856
+
+
+plot(IC7)
+
+system.time(IC8 <- radiusMinimaxIC(L2Fam=N, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0.01, upRad=1.8))
+# user system elapsed
+# 565.58 0.21 586.05
+
+IC8
+
+#An object of class “TotalVarIC”
+#### name: IC of total variation type
+#
+#### L2-differentiable parametric family: Negative 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.2944932
+#
+#### clipLo: [1] -0.06497353
+#### clipUp: [1] 0.05431373
+#### stand:
+# prob
+#prob 0.003432502
+#
+#### Infos:
+# method message
+#[1,] "radiusMinimaxIC" "radius minimax IC for radius interval [0.01, 1.8]"
+#[2,] "radiusMinimaxIC" "least favorable radius: 0.294"
+#[3,] "radiusMinimaxIC" "maximum ‘asMSE’-inefficiency: 1.168"
+
+
+checkIC(IC8)
+
+#precision of centering: 6.400051e-12
+#precision of Fisher consistency:
+# prob
+#prob -1.404645e-09
+#maximum deviation
+# 1.404645e-09
+
+
+Risks(IC8)
+
+#$asCov
+# prob
+#prob 0.003987582
+#
+#$asBias
+#$asBias$value
+#[1] 0.1192873
+#
+#$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.003987582
+#
+#$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.005221649
+#
+#$asMSE$r
+#[1] 0.2944932
+#
+#$asMSE$at
+#An object of class “TotalVarNeighborhood”
+#type: (uncond.) total variation neighborhood
+#radius: 0.2944932
+
+plot(IC8)
+
+
+#-------------------------------------------------------------------------------
+## least favorable radius
+#-------------------------------------------------------------------------------
+## (may take quite some time!)
+system.time(r.rho1 <- leastFavorableRadius(L2Fam=N, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5))
+
+#current radius: 0.3820278 inefficiency: 1.040429
+#current radius: 0.6180722 inefficiency: 1.044464
+#current radius: 0.7639556 inefficiency: 1.041844
+#current radius: 0.5931816 inefficiency: 1.044631
+#current radius: 0.5546088 inefficiency: 1.044691
+#current radius: 0.5644109 inefficiency: 1.044698
+#current radius: 0.5640279 inefficiency: 1.0447
+#current radius: 0.5599945 inefficiency: 1.044697
+#current radius: 0.5624873 inefficiency: 1.044701
+#current radius: 0.5631909 inefficiency: 1.044701
+#current radius: 0.5627771 inefficiency: 1.044701
+#current radius: 0.5625595 inefficiency: 1.044701
+#current radius: 0.5626002 inefficiency: 1.044701
+#current radius: 0.5625595 inefficiency: 1.044701
+# user system elapsed
+# 361.25 0.12 369.05
+
+## same as for binomial????
+
+r.rho1
+
+#$rho
+#[1] 0.5
+#
+#$leastFavorableRadius
+#[1] 0.5625595
+#
+#$`asMSE-inefficiency`
+#[1] 1.044701
+
+system.time(r.rho2 <- leastFavorableRadius(L2Fam=N, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=0.5))
+
+#current radius: 0.3820278 inefficiency: 1.041727
+#current radius: 0.6180722 inefficiency: 1.027733
+#current radius: 0.2361444 inefficiency: 1.043317
+#current radius: 0.2660735 inefficiency: 1.044275
+#current radius: 0.2943633 inefficiency: 1.044409
+#current radius: 0.2852759 inefficiency: 1.04443
+#current radius: 0.2866889 inefficiency: 1.044456
+#current radius: 0.2893884 inefficiency: 1.044426
+#current radius: 0.2872589 inefficiency: 1.044427
+#current radius: 0.2862418 inefficiency: 1.044439
+#current radius: 0.2869066 inefficiency: 1.044442
+#current radius: 0.2865879 inefficiency: 1.044429
+#current radius: 0.2867296 inefficiency: 1.044453
+#current radius: 0.2866482 inefficiency: 1.044425
+#current radius: 0.2866889 inefficiency: 1.044456
+# user system elapsed
+# 4891.07 1.90 5063.44
+
+r.rho2
+
+#$rho
+#[1] 0.5
+#
+#$leastFavorableRadius
+#[1] 0.2866889
+#
+#$`asMSE-inefficiency`
+#[1] 1.044456
+
+
+###############################################################################
+## k-step (k >= 1) estimation
+################################################################################
+
+## one-step estimation
+## 1. generate a contaminated sample
+ind <- rbinom(100, size=1, prob=0.05)
+x <- rnbinom(100, size=25, prob=(1-ind)*0.25 + ind*0.01)
+
+### MLE:
+
+(estML <- MLEstimator(x=x, NbinomFamily(size=25)))
+
+#Evaluations of Maximum likelihood estimate:
+#-------------------------------------------
+#An object of class “Estimate”
+#generated by call
+# MLEstimator(x = x, ParamFamily = NbinomFamily(size = 25))
+#samplesize: 100
+#estimate:
+# prob
+# 0.135624871
+# (0.002521857)
+#fixed part of the parameter:
+#size
+# 25
+#asymptotic (co)variance (multiplied with samplesize):
+#[1] 0.0006359763
+#Criterion:
+#negative log-likelihood
+# 1868.396
+
+## 2. Kolmogorov(-Smirnov) minimum distance estimator
+
+(est0 <- MDEstimator(x=x, NbinomFamily(size=25)))
+
+#Evaluations of Minimum Kolmogorov distance estimate:
+#----------------------------------------------------
+#An object of class “Estimate”
+#generated by call
+# MDEstimator(x = x, ParamFamily = NbinomFamily(size = 25))
+#samplesize: 100
+#estimate:
+# prob
+#0.2471440
+#fixed part of the parameter:
+#size
+# 25
+#Criterion:
+#Kolmogorov distance
+# 0.05226461
+
+### 3.1. one-step estimation: radius known
+
+### ac) Define infinitesimal robust model
+RobN3 <- InfRobModel(center=NbinomFamily(size=25, prob=estimate(est0)),
+ neighbor=ContNeighborhood(radius=0.5))
+
+## bc) Compute optimally robust IC
+
+IC9 <- optIC(model=RobN3, risk=asMSE())
+checkIC(IC9)
+
+#precision of centering: 5.93858e-07
+#precision of Fisher consistency:
+# prob
+#prob 8.014067e-05
+#maximum deviation
+# 8.014067e-05
+
+(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.251879937
+# (0.004632576)
+#fixed part of the parameter:
+#size
+# 25
+#asymptotic (co)variance (multiplied with samplesize):
+#[1] 0.002146076
+#Infos:
+# method
+#[1,] "oneStepEstimator"
+#[2,] "oneStepEstimator"
+# message
+#[1,] "1-step estimate for Negative Binomial family"
+#[2,] "computation of IC, trafo, asvar and asbias via useLast = FALSE"
+#asymptotic bias:
+#[1] 0.03064529
+#(partial) influence curve:
+#An object of class “ContIC”
+#### name: IC of contamination type
+#
+#### L2-differentiable parametric family: Negative Binomial family
+#### param: An object of class "ParamFamParameter"
+#name: probability of success
+#prob: 0.249199096954169
+#fixed part of param.:
+# size: 25
+#trafo:
+# prob
+#prob 1
+#
+#### neighborhood radius: 0.5
+#
+#### clip: [1] 0.06129058
+#### cent: [1] 0.003713275
+#### stand:
+# prob
+#prob 0.00308521
+#
+#### Infos:
+# method message
+#[1,] "optIC" "optimally robust IC for ‘asMSE’"
+#steps:
+#[1] 1
+
+est1c1 <- roptest(x, NbinomFamily(size = 25), eps = 0.05, initial.est = est0)
+checkIC(pIC(est1c1))
+
+#precision of centering: 3.376922e-18
+#precision of Fisher consistency:
+# prob
+#prob 7.04511e-05
+#maximum deviation
+# 7.04511e-05
+
+est1c2 <- roptest(x, NbinomFamily(size = 25), eps = 0.05, distance = KolmogorovDist)
+checkIC(pIC(est1c2))
+
+#precision of centering: 3.376922e-18
+#precision of Fisher consistency:
+# prob
+#prob 7.04511e-05
+#maximum deviation
+# 7.04511e-05
+
+est1c3 <- roptest(x, NbinomFamily(size = 25), eps = 0.025)
+checkIC(pIC(est1c3))
+
+#precision of centering: 8.748485e-11
+#precision of Fisher consistency:
+# prob
+#prob 8.488734e-05
+#maximum deviation
+# 8.488734e-05
+
+
+estimate(est1c)
+# prob
+#0.2518799
+estimate(est1c1)
+# prob
+#0.2518792
+estimate(est1c2)
+# prob
+#0.2518792
+estimate(est1c3)
+# prob
+#0.2502157
+confint(est1c, symmetricBias())
+#A[n] asymptotic (LAN-based), uniform (bias-aware)
+# confidence interval:
+#for symmetric Bias
+# 2.5 % 97.5 %
+#prob 0.2426583 0.2611016
+#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 407
More information about the Robast-commits
mailing list