[Robast-commits] r46 - in pkg/ROptEst/inst: . scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 19 07:09:45 CET 2008
Author: stamats
Date: 2008-02-19 07:09:45 +0100 (Tue, 19 Feb 2008)
New Revision: 46
Added:
pkg/ROptEst/inst/scripts/
pkg/ROptEst/inst/scripts/BinomialModel.R
pkg/ROptEst/inst/scripts/ExponentialScaleModel.R
pkg/ROptEst/inst/scripts/GammaModel.R
pkg/ROptEst/inst/scripts/GumbelLocationModel.R
pkg/ROptEst/inst/scripts/LognormalAndNormalModel.R
pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
pkg/ROptEst/inst/scripts/NormalScaleModel.R
pkg/ROptEst/inst/scripts/PoissonModel.R
pkg/ROptEst/inst/scripts/UnderOverShootRisk.R
Log:
added folder scripts again
Added: pkg/ROptEst/inst/scripts/BinomialModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/BinomialModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/BinomialModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,132 @@
+###############################################################################
+## Example: Binomial Family
+###############################################################################
+require(ROptEst)
+
+## generates Binomial Family with
+## m = 25 and probability of success theta = 0.25
+B <- BinomFamily(size = 25, prob = 0.25)
+B # show B
+plot(B) # plot of Binom(size = 25, prob = 0.25) and L_2 derivative
+checkL2deriv(B)
+
+## classical optimal IC
+IC0 <- optIC(model = B, risk = asCov())
+IC0 # show IC
+plot(IC0) # plot IC
+checkIC(IC0)
+Risks(IC0)
+
+## lower case radius
+lowerCaseRadius(L2Fam = B, neighbor = ContNeighborhood(), risk = asMSE())
+lowerCaseRadius(L2Fam = B, neighbor = TotalVarNeighborhood(), risk = asMSE())
+
+## L_2 family + infinitesimal neighborhood
+RobB1 <- InfRobModel(center = B, neighbor = ContNeighborhood(radius = 0.5))
+RobB1 # show RobB1
+(RobB2 <- InfRobModel(center = B, neighbor = TotalVarNeighborhood(radius = 0.5)))
+
+## MSE solution
+system.time(IC1 <- optIC(model=RobB1, risk=asMSE()), gcFirst = TRUE)
+IC1
+checkIC(IC1)
+Risks(IC1)
+getRiskIC(IC1, asBias(), ContNeighborhood()) # standardized bias
+getRiskIC(IC1, asMSE(), ContNeighborhood(radius = 0.5))
+
+(Cov1 <- getRiskIC(IC1, asCov()))
+(mse1 <- getRiskIC(IC1, asMSE(), TotalVarNeighborhood(radius = 0.5)))
+(bias1 <- getRiskIC(IC1, asBias(), TotalVarNeighborhood()))
+## only suboptimal -> ToDo-List
+addRisk(IC1) <- list(Cov1, mse1, bias1)
+Risks(IC1)
+plot(IC1)
+
+system.time(IC2 <- optIC(model=RobB2, risk=asMSE()), gcFirst = TRUE)
+IC2
+checkIC(IC2)
+Risks(IC2)
+getRiskIC(IC2, asMSE(), TotalVarNeighborhood(radius = 0.5))
+getRiskIC(IC2, asBias(), TotalVarNeighborhood())
+getRiskIC(IC2, asBias(), ContNeighborhood())
+Cov2 <- getRiskIC(IC2, asCov())
+addRisk(IC2) <- Cov2
+Risks(IC2)
+plot(IC2)
+
+## lower case solutions
+(IC3 <- optIC(model=RobB1, risk=asBias()))
+checkIC(IC3)
+Risks(IC3)
+plot(IC3)
+
+(IC4 <- optIC(model=RobB2, risk=asBias()))
+checkIC(IC4)
+Risks(IC4)
+plot(IC4)
+
+
+## Hampel solution
+(IC5 <- optIC(model=RobB1, risk=asHampel(bound=clip(IC1))))
+checkIC(IC5)
+Risks(IC5)
+plot(IC5)
+
+(IC6 <- optIC(model=RobB2, risk=asHampel(bound=Risks(IC2)$asBias), maxiter = 200))
+checkIC(IC6)
+Risks(IC6)
+plot(IC6)
+
+
+## radius minimax IC
+system.time(IC7 <- radiusMinimaxIC(L2Fam=B, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=1), gcFirst = TRUE)
+IC7
+checkIC(IC7)
+Risks(IC7)
+plot(IC7)
+
+system.time(IC8 <- radiusMinimaxIC(L2Fam=B, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=1), gcFirst = TRUE)
+IC8
+checkIC(IC8)
+Risks(IC8)
+plot(IC8)
+
+
+## least favorable radius
+## (may take quite some time!)
+system.time(r.rho1 <- leastFavorableRadius(L2Fam=B, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5), gcFirst = TRUE)
+r.rho1
+system.time(r.rho2 <- leastFavorableRadius(L2Fam=B, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=0.5), gcFirst = TRUE)
+r.rho2
+
+## 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 <- ksEstimator(x=x, Binom(size=25), param = "prob"))
+
+## 3. one-step estimation: radius known
+RobB3 <- InfRobModel(center=BinomFamily(size=25, prob=est0$prob),
+ neighbor=ContNeighborhood(radius=0.5))
+IC9 <- optIC(model=RobB3, risk=asMSE())
+(est1 <- oneStepEstimator(x, IC=IC9, start=est0$prob))
+
+RobB4 <- InfRobModel(center=BinomFamily(size=25, prob=est0$prob),
+ neighbor=TotalVarNeighborhood(radius=0.25))
+IC10 <- optIC(model=RobB4, risk=asMSE())
+(est1 <- oneStepEstimator(x, IC=IC10, start=est0$prob))
+
+## 4. one-step estimation: radius interval
+IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=est0$prob),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(est2 <- oneStepEstimator(x, IC=IC11, start=est0$prob))
+
+IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=est0$prob),
+ neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(est2 <- oneStepEstimator(x, IC=IC12, start=est0$prob))
Added: pkg/ROptEst/inst/scripts/ExponentialScaleModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/ExponentialScaleModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/ExponentialScaleModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,88 @@
+###############################################################################
+## Example: Exponential Scale Family
+###############################################################################
+require(ROptEst)
+
+## generates Exponential Scale Family with rate = 1
+E1 <- ExpScaleFamily(rate = 1)
+E1 # show E1
+plot(E1) # plot of Exp(rate = 1) and L_2 derivative
+checkL2deriv(E1)
+
+# classical optimal IC
+E1.IC0 <- optIC(model = E1, risk = asCov())
+E1.IC0 # show IC
+checkIC(E1.IC0)
+Risks(E1.IC0)
+plot(E1.IC0) # plot IC
+
+# L_2 family + infinitesimal neighborhood
+E1.Rob1 <- InfRobModel(center = E1, neighbor = ContNeighborhood(radius = 0.5))
+E1.Rob1 # show E1.Rob1
+E1.Rob2 <- InfRobModel(center = E1, neighbor = TotalVarNeighborhood(radius = 0.5))
+
+# MSE solution
+(E1.IC1 <- optIC(model=E1.Rob1, risk=asMSE()))
+checkIC(E1.IC1)
+Risks(E1.IC1)
+plot(E1.IC1)
+(E1.IC2 <- optIC(model=E1.Rob2, risk=asMSE()))
+checkIC(E1.IC2)
+Risks(E1.IC2)
+plot(E1.IC2)
+
+# lower case solutions
+(E1.IC3 <- optIC(model=E1.Rob1, risk=asBias()))
+checkIC(E1.IC3)
+Risks(E1.IC3)
+plot(E1.IC3)
+(E1.IC4 <- optIC(model=E1.Rob2, risk=asBias()))
+checkIC(E1.IC4)
+Risks(E1.IC4)
+plot(E1.IC4)
+
+# Hampel solution
+(E1.IC5 <- optIC(model=E1.Rob1, risk=asHampel(bound=clip(E1.IC1))))
+checkIC(E1.IC5)
+Risks(E1.IC5)
+plot(E1.IC5)
+(E1.IC6 <- optIC(model=E1.Rob2, risk=asHampel(bound=Risks(E1.IC2)$asBias), maxiter = 200))
+checkIC(E1.IC6)
+Risks(E1.IC6)
+plot(E1.IC6)
+
+# radius minimax IC
+(E1.IC7 <- radiusMinimaxIC(L2Fam=E1, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=0.5))
+checkIC(E1.IC7)
+Risks(E1.IC7)
+(E1.IC8 <- radiusMinimaxIC(L2Fam=E1, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=0.5))
+checkIC(E1.IC8)
+Risks(E1.IC8)
+
+# least favorable radius
+# (may take quite some time!)
+(E1.r.rho1 <- leastFavorableRadius(L2Fam=E1, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5))
+(E1.r.rho2 <- leastFavorableRadius(L2Fam=E1, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=1/3))
+
+## one-step estimation
+## 1. generate a contaminated sample
+ind <- rbinom(1e2, size=1, prob=0.05)
+E1.x <- rexp(1e2, rate=(1-ind)*2+ind*10)
+
+## 2. Kolmogorov(-Smirnov) minimum distance estimator
+(E1.est0 <- ksEstimator(x=E1.x, Exp()))
+
+## 3. one-step estimation: radius known
+E1.Rob3 <- InfRobModel(center=ExpScaleFamily(rate=E1.est0$rate),
+ neighbor=ContNeighborhood(radius=0.5))
+E1.IC9 <- optIC(model=E1.Rob3, risk=asMSE())
+(E1.est1 <- oneStepEstimator(E1.x, IC=E1.IC9, start=E1.est0$rate))
+
+## 4. one-step estimation: radius interval
+E1.IC10 <- radiusMinimaxIC(L2Fam=ExpScaleFamily(rate=E1.est0$rate),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(E1.est2 <- oneStepEstimator(E1.x, IC=E1.IC10, start=E1.est0$rate))
Added: pkg/ROptEst/inst/scripts/GammaModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/GammaModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/GammaModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,73 @@
+###############################################################################
+## Example: Gamma Family
+###############################################################################
+require(ROptEst)
+
+## generates Gamma Family with
+## scale = 1 and shape = 1
+G <- GammaFamily(scale = 1, shape = 2)
+G # show G
+plot(G) # plot of Gammad(scale = 1, shape = 2) and L_2 derivative
+checkL2deriv(G)
+
+## classical optimal IC
+IC0 <- optIC(model = G, risk = asCov())
+IC0 # show IC
+system.time(checkIC(IC0), gcFirst = TRUE)
+Risks(IC0)
+plot(IC0) # plot IC
+
+## L_2 family + infinitesimal neighborhood
+RobG1 <- InfRobModel(center = G, neighbor = ContNeighborhood(radius = 0.5))
+RobG1 # show RobB1
+
+## MSE solution
+system.time(IC1 <- optIC(model=RobG1, risk=asMSE()), gcFirst = TRUE)
+IC1
+checkIC(IC1)
+Risks(IC1)
+plot(IC1)
+x11()
+infoPlot(IC1)
+
+## lower case solutions
+system.time(IC2 <- optIC(model=RobG1, risk=asBias(), tol = 1e-10), gcFirst = TRUE)
+IC2
+checkIC(IC2)
+Risks(IC2)
+plot(IC2)
+x11()
+infoPlot(IC2)
+
+## Hampel solution
+system.time(IC3 <- optIC(model=RobG1, risk=asHampel(bound=clip(IC1))), gcFirst = TRUE)
+IC3
+checkIC(IC3)
+Risks(IC3)
+plot(IC3)
+x11()
+infoPlot(IC3)
+
+## radius minimax IC
+## takes quite some time - about 30 min.
+system.time(IC4 <- radiusMinimaxIC(L2Fam=G, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=Inf), gcFirst = TRUE)
+
+## least favorable radius
+## takes quite some time - several hours!
+#system.time(r.rho1 <- leastFavorableRadius(L2Fam=G, neighbor=ContNeighborhood(),
+# risk=asMSE(), rho=0.5))
+
+## one-step estimation
+## 1. generate a contaminated sample
+ind <- rbinom(100, size=1, prob=0.05)
+x <- (1-ind)*rgamma(100, scale = 1, shape = 2) + ind*10
+
+## 2. Kolmogorov(-Smirnov) minimum distance estimator
+(est0 <- ksEstimator(x=x, Gammad()))
+
+## 3. one-step estimation: radius known
+RobG3 <- InfRobModel(center=GammaFamily(scale = est0$scale, shape = est0$shape),
+ neighbor=ContNeighborhood(radius=0.5))
+IC9 <- optIC(model=RobG3, risk=asMSE())
+(est1 <- oneStepEstimator(x, IC=IC9, start=est0))
Added: pkg/ROptEst/inst/scripts/GumbelLocationModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/GumbelLocationModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/GumbelLocationModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,135 @@
+###############################################################################
+## Example: Gumbel Location Family
+## computations numerically less stable than in case of the
+## Exponential Scale Family
+###############################################################################
+require(ROptEst)
+
+## generates Gumbel Location Family with loc = 0
+## (known scale = 1)
+distrExOptions(ElowerTruncQuantile, 1e-15) # non-finite function value in integrate
+G0 <- GumbelLocationFamily(loc=0, scale=1)
+G0 # show G0
+plot(G0) # plot of Gumbel(loc = 0, scale = 1) and L_2 derivative
+checkL2deriv(G0)
+
+# classical optimal IC
+G0.IC0 <- optIC(model = G0, risk = asCov())
+G0.IC0 # show IC
+plot(G0.IC0) # plot IC
+checkIC(G0.IC0)
+Risks(G0.IC0)
+
+# L_2 family + infinitesimal neighborhood
+G0.Rob1 <- InfRobModel(center = G0, neighbor = ContNeighborhood(radius = 0.5))
+G0.Rob1 # show G0.Rob1
+G0.Rob2 <- InfRobModel(center = G0, neighbor = TotalVarNeighborhood(radius = 0.5))
+
+# MSE solution
+E1.Rob1 <- InfRobModel(center = ExpScaleFamily(), neighbor = ContNeighborhood(radius = 0.5))
+(E1.IC1 <- optIC(model=E1.Rob1, risk=asMSE()))
+G0.IC1 <- optIC(model=G0.Rob1, risk=asMSE())
+checkIC(G0.IC1)
+Risks(G0.IC1)
+clip(E1.IC1)
+cent(E1.IC1)
+stand(E1.IC1)
+clip(G0.IC1)
+cent(G0.IC1)
+stand(G0.IC1)
+# alternatively
+G0.IC11 <- E1.IC1 # rate = 1!
+CallL2Fam(G0.IC11) <- call("GumbelLocationFamily")
+cent(G0.IC11) <- -cent(E1.IC1)
+G0.IC11
+checkIC(G0.IC11)
+Risks(G0.IC11)
+
+E1.Rob2 <- InfRobModel(center = ExpScaleFamily(), neighbor = TotalVarNeighborhood(radius = 0.5))
+E1.IC2 <- optIC(model=E1.Rob2, risk=asMSE())
+#distrExOptions(ElowerTruncQuantile, 1e-15)
+G0.IC2 <- optIC(model=G0.Rob2, risk=asMSE())
+checkIC(G0.IC2)
+Risks(G0.IC2)
+clipLo(E1.IC2)
+clipUp(E1.IC2)
+stand(E1.IC2)
+clipLo(G0.IC2)
+clipUp(G0.IC2)
+stand(G0.IC2)
+# alternatively
+G0.IC21 <- E1.IC2 # rate = 1!
+CallL2Fam(G0.IC21) <- call("GumbelLocationFamily")
+clipLo(G0.IC21) <- -clipUp(E1.IC2)
+clipUp(G0.IC21) <- -clipLo(E1.IC2)
+G0.IC21
+checkIC(G0.IC21)
+Risks(G0.IC21)
+
+# lower case solutions
+(G0.IC3 <- optIC(model=G0.Rob1, risk=asBias()))
+checkIC(G0.IC3)
+Risks(G0.IC3)
+(G0.IC4 <- optIC(model=G0.Rob2, risk=asBias()))
+checkIC(G0.IC4)
+Risks(G0.IC4)
+
+# Hampel solution
+(G0.IC5 <- optIC(model=G0.Rob1, risk=asHampel(bound=clip(G0.IC1))))
+checkIC(G0.IC5)
+Risks(G0.IC5)
+(G0.IC6 <- optIC(model=G0.Rob2, risk=asHampel(bound=Risks(G0.IC2)$asBias), maxiter = 100))
+checkIC(G0.IC6)
+Risks(G0.IC6)
+
+# radius minimax IC
+# numerically instable for small 'loRad'!
+# => use connection to ExpScaleFamily for computations
+#(G0.IC7 <- radiusMinimaxIC(L2Fam=G0, neighbor=ContNeighborhood(),
+# risk=asMSE(), loRad=0.5, upRad=1.0))
+#checkIC(G0.IC7)
+#Risks(G0.IC7)
+#(G0.IC8 <- radiusMinimaxIC(L2Fam=G0, neighbor=TotalVarNeighborhood(),
+# risk=asMSE(), loRad=0.5, upRad=1.0))
+#checkIC(G0.IC8)
+#Risks(G0.IC8)
+
+# least favorable radius
+# numerically instable!
+# => use connection to ExpScaleFamily for computations
+#(G0.r.rho1 <- leastFavorableRadius(L2Fam=G0, neighbor=ContNeighborhood(),
+# risk=asMSE(), rho=0.5))
+#(G0.r.rho2 <- leastFavorableRadius(L2Fam=G0, neighbor=TotalVarNeighborhood(),
+# risk=asMSE(), rho=1/3))
+
+## one-step estimation
+## 1. generate a contaminated sample
+ind <- rbinom(1e2, size=1, prob=0.05)
+G0.x <- rgumbel(1e2, loc=(1-ind)*0.5+ind*1)
+
+## 2. Kolmogorov(-Smirnov) minimum distance estimator
+(G0.est0 <- ksEstimator(x=G0.x, Gumbel(), param = "loc"))
+
+## 3. one-step estimation: radius known
+G0.Rob3 <- InfRobModel(center=GumbelLocationFamily(loc=G0.est0$loc),
+ neighbor=ContNeighborhood(radius=0.5))
+G0.IC9 <- optIC(model=G0.Rob3, risk=asMSE())
+(G0.est1 <- oneStepEstimator(G0.x, IC=G0.IC9, start=G0.est0$loc))
+
+## 4. M estimation: radius known
+G0.Rob31 <- InfRobModel(center=GumbelLocationFamily(loc=0),
+ neighbor=ContNeighborhood(radius=0.5))
+G0.IC91 <- optIC(model=G0.Rob31, risk=asMSE())
+(G0.est11 <- locMEstimator(G0.x, IC=G0.IC91))
+
+## 5. one-step estimation: radius interval
+G0.IC10 <- radiusMinimaxIC(L2Fam=GumbelLocationFamily(loc=G0.est0$loc),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0.5, upRad=1)
+(G0.est2 <- oneStepEstimator(G0.x, IC=G0.IC10, start=G0.est0$loc))
+
+## 6. M estimation: radius interval
+G0.IC101 <- radiusMinimaxIC(L2Fam=GumbelLocationFamily(),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0.5, upRad=1)
+(G0.est21 <- locMEstimator(G0.x, IC=G0.IC101))
+
+distrExOptions(ElowerTruncQuantile, 0) # default
Added: pkg/ROptEst/inst/scripts/LognormalAndNormalModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/LognormalAndNormalModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/LognormalAndNormalModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,155 @@
+###############################################################################
+## Example: Lognormal Scale and Normal Location
+###############################################################################
+require(ROptEst)
+
+## generates Lognormal Scale Family with rate = 1
+LN1 <- LnormScaleFamily()
+LN1 # show LN1
+plot(LN1) # plot of Exp(rate = 1) and L_2 derivative
+checkL2deriv(LN1)
+
+## generates Normal Location Family with mean = 0
+N0 <- NormLocationFamily(mean=0, sd=1)
+N0 # show G0
+plot(N0) # plot of Norm(mean = 0, sd = 1) and L_2 derivative
+checkL2deriv(N0)
+
+
+# classical optimal IC
+LN1.IC0 <- optIC(model = LN1, risk = asCov())
+LN1.IC0 # show IC
+plot(LN1.IC0) # plot IC
+checkIC(LN1.IC0)
+Risks(LN1.IC0)
+N0.IC0 <- optIC(model = N0, risk = asCov())
+N0.IC0 # show IC
+plot(N0.IC0) # plot IC
+checkIC(N0.IC0)
+Risks(N0.IC0)
+
+
+# L_2 family + infinitesimal neighborhood
+LN1.Rob1 <- InfRobModel(center = LN1, neighbor = ContNeighborhood(radius = 0.5))
+LN1.Rob1 # show LN1.Rob1
+LN1.Rob2 <- InfRobModel(center = LN1, neighbor = TotalVarNeighborhood(radius = 0.25))
+N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 0.5))
+N0.Rob1 # show N0.Rob1
+N0.Rob2 <- InfRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = 0.25))
+
+
+# MSE solution
+LN1.IC1 <- optIC(model=LN1.Rob1, risk=asMSE())
+checkIC(LN1.IC1)
+Risks(LN1.IC1)
+plot(LN1.IC1)
+
+N0.IC1 <- optIC(model=N0.Rob1, risk=asMSE())
+checkIC(N0.IC1)
+Risks(N0.IC1)
+plot(N0.IC1)
+
+clip(LN1.IC1)
+cent(LN1.IC1)
+stand(LN1.IC1)
+clip(N0.IC1)
+cent(N0.IC1)
+stand(N0.IC1)
+
+LN1.IC2 <- optIC(model=LN1.Rob2, risk=asMSE())
+checkIC(LN1.IC2)
+Risks(LN1.IC2)
+plot(LN1.IC2)
+
+N0.IC2 <- optIC(model=N0.Rob2, risk=asMSE())
+checkIC(N0.IC2)
+Risks(N0.IC2)
+plot(N0.IC2)
+
+clipLo(LN1.IC2)
+clipUp(LN1.IC2)
+stand(LN1.IC2)
+clipLo(N0.IC2)
+clipUp(N0.IC2)
+stand(N0.IC2)
+
+
+# lower case solutions
+LN1.IC3 <- optIC(model=LN1.Rob1, risk=asBias())
+checkIC(LN1.IC3)
+Risks(LN1.IC3)
+plot(LN1.IC3)
+
+N0.IC3 <- optIC(model=N0.Rob1, risk=asBias())
+checkIC(N0.IC3)
+Risks(N0.IC3)
+plot(N0.IC3)
+
+LN1.IC4 <- optIC(model=LN1.Rob2, risk=asBias())
+checkIC(LN1.IC4)
+Risks(LN1.IC4)
+plot(LN1.IC4)
+
+N0.IC4 <- optIC(model=N0.Rob2, risk=asBias())
+checkIC(N0.IC4)
+Risks(N0.IC4)
+plot(N0.IC4)
+
+
+# Hampel solution
+LN1.IC5 <- optIC(model=LN1.Rob1, risk=asHampel(bound=clip(LN1.IC1)))
+checkIC(LN1.IC5)
+Risks(LN1.IC5)
+plot(LN1.IC5)
+
+N0.IC5 <- optIC(model=N0.Rob1, risk=asHampel(bound=clip(N0.IC1)))
+checkIC(N0.IC5)
+Risks(N0.IC5)
+plot(N0.IC5)
+
+LN1.IC6 <- optIC(model=LN1.Rob2, risk=asHampel(bound=Risks(LN1.IC2)$asBias))
+checkIC(LN1.IC6)
+Risks(LN1.IC6)
+plot(LN1.IC6)
+
+N0.IC6 <- optIC(model=N0.Rob2, risk=asHampel(bound=Risks(N0.IC2)$asBias))
+checkIC(N0.IC6)
+Risks(N0.IC6)
+plot(N0.IC6)
+
+# radius minimax IC
+(LN1.IC7 <- radiusMinimaxIC(L2Fam=LN1, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=0.5))
+checkIC(LN1.IC7)
+Risks(LN1.IC7)
+plot(LN1.IC7)
+
+(N0.IC7 <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0.1, upRad=0.5))
+checkIC(N0.IC7)
+Risks(N0.IC7)
+plot(N0.IC7)
+
+(LN1.IC8 <- radiusMinimaxIC(L2Fam=LN1, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=0.25))
+checkIC(LN1.IC8)
+Risks(LN1.IC8)
+plot(LN1.IC8)
+
+(N0.IC8 <- radiusMinimaxIC(L2Fam=N0, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=0.25))
+checkIC(N0.IC8)
+Risks(N0.IC8)
+plot(N0.IC8)
+
+
+# least favorable radius
+# (may take quite some time!)
+(LN1.r.rho1 <- leastFavorableRadius(L2Fam=LN1, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5))
+(N0.r.rho1 <- leastFavorableRadius(L2Fam=N0, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5))
+(LN1.r.rho2 <- leastFavorableRadius(L2Fam=LN1, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=1/3))
+(N0.r.rho2 <- leastFavorableRadius(L2Fam=N0, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=1/3))
Added: pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,79 @@
+###############################################################################
+## Example: Normal location and scale
+###############################################################################
+require(ROptEst)
+
+## generates normal location and scale family with mean = 0 and sd = 1
+N0 <- NormLocationScaleFamily(mean=0, sd=1)
+N0 # show G0
+plot(N0) # plot of Norm(mean = 0, sd = 1) and L_2 derivative
+checkL2deriv(N0)
+
+# classical optimal IC
+N0.IC0 <- optIC(model = N0, risk = asCov())
+N0.IC0 # show IC
+checkIC(N0.IC0)
+Risks(N0.IC0)
+plot(N0.IC0) # plot IC
+
+# L_2 family + infinitesimal neighborhood
+N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 0.5))
+N0.Rob1 # show N0.Rob1
+
+# MSE solution
+system.time(N0.IC1 <- optIC(model = N0.Rob1, risk = asMSE()), gcFirst = TRUE)
+checkIC(N0.IC1)
+Risks(N0.IC1)
+plot(N0.IC1)
+infoPlot(N0.IC1)
+
+# lower case solutions
+(N0.IC2 <- optIC(model = N0.Rob1, risk = asBias(), tol = 1e-10))
+checkIC(N0.IC2)
+Risks(N0.IC2)
+plot(N0.IC2)
+infoPlot(N0.IC2)
+
+# Hampel solution
+(N0.IC3 <- optIC(model = N0.Rob1, risk = asHampel(bound = clip(N0.IC1))))
+checkIC(N0.IC3)
+Risks(N0.IC3)
+plot(N0.IC3)
+infoPlot(N0.IC3)
+
+# radius minimax IC
+# (may take quite some time!)
+(N0.IC4 <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=Inf))
+checkIC(N0.IC4)
+Risks(N0.IC4)
+plot(N0.IC4)
+infoPlot(N0.IC4)
+
+# least favorable radius
+# (may take quite some time!)
+#N0.r.rho1 <- leastFavorableRadius(L2Fam=N0, neighbor=ContNeighborhood(),
+# risk=asMSE(), rho=0.5)
+
+## one-step estimation
+## 1. generate a contaminated sample
+ind <- rbinom(100, size=1, prob=0.05)
+x <- rnorm(100, mean=0, sd=(1-ind) + ind*9)
+
+## 2. Kolmogorov(-Smirnov) minimum distance estimator
+(est0 <- ksEstimator(x=x, Norm()))
+
+## 3. one-step estimation: radius known
+N1 <- NormLocationScaleFamily(mean=est0$mean, sd=est0$sd)
+N1.Rob <- InfRobModel(center = N1, neighbor = ContNeighborhood(radius = 0.5))
+IC1 <- optIC(model = N1.Rob, risk = asMSE())
+(est1 <- oneStepEstimator(x, IC1, est0))
+
+## 4. one-step estimation: radius unknown
+## rough estimate: 1-10% contamination
+## => r\in[0.1,1.0]
+
+## takes some time
+IC2 <- radiusMinimaxIC(L2Fam=N1, neighbor=ContNeighborhood(),risk=asMSE(),
+ loRad=0.1, upRad=1.0)
+(est2 <- oneStepEstimator(x, IC2, est0))
Added: pkg/ROptEst/inst/scripts/NormalScaleModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/NormalScaleModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/NormalScaleModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,72 @@
+###############################################################################
+## Example: Normal Scale
+###############################################################################
+require(ROptEst)
+
+## generates Normal Scale Family with scale = 1
+N0 <- NormScaleFamily(mean=0, sd=1)
+N0 # show G0
+plot(N0) # plot of Norm(mean = 0, sd = 1) and L_2 derivative
+checkL2deriv(N0)
+
+# classical optimal IC
+N0.IC0 <- optIC(model = N0, risk = asCov())
+N0.IC0 # show IC
+plot(N0.IC0) # plot IC
+checkIC(N0.IC0)
+Risks(N0.IC0)
+
+# L_2 family + infinitesimal neighborhood
+N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 0.5))
+N0.Rob1 # show N0.Rob1
+N0.Rob2 <- InfRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = 0.5))
+
+# MSE solution
+(N0.IC1 <- optIC(model=N0.Rob1, risk=asMSE()))
+checkIC(N0.IC1)
+Risks(N0.IC1)
+plot(N0.IC1)
+
+(N0.IC2 <- optIC(model=N0.Rob2, risk=asMSE()))
+checkIC(N0.IC2)
+Risks(N0.IC2)
+plot(N0.IC2)
+
+# lower case solutions
+(N0.IC3 <- optIC(model=N0.Rob1, risk=asBias()))
+checkIC(N0.IC3)
+Risks(N0.IC3)
+plot(N0.IC3)
+(N0.IC4 <- optIC(model=N0.Rob2, risk=asBias()))
+checkIC(N0.IC4)
+Risks(N0.IC4)
+plot(N0.IC4)
+
+# Hampel solution
+(N0.IC5 <- optIC(model=N0.Rob1, risk=asHampel(bound=clip(N0.IC1))))
+checkIC(N0.IC5)
+Risks(N0.IC5)
+plot(N0.IC5)
+(N0.IC6 <- optIC(model=N0.Rob2, risk=asHampel(bound=Risks(N0.IC2)$asBias), maxiter = 200))
+checkIC(N0.IC6)
+Risks(N0.IC6)
+plot(N0.IC6)
+
+# radius minimax IC
+(N0.IC7 <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=Inf))
+checkIC(N0.IC7)
+Risks(N0.IC7)
+plot(N0.IC7)
+(N0.IC8 <- radiusMinimaxIC(L2Fam=N0, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=Inf))
+checkIC(N0.IC8)
+Risks(N0.IC8)
+plot(N0.IC8)
+
+# least favorable radius
+# (may take quite some time!)
+(N0.r.rho1 <- leastFavorableRadius(L2Fam=N0, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5))
+(N0.r.rho2 <- leastFavorableRadius(L2Fam=N0, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=1/3))
Added: pkg/ROptEst/inst/scripts/PoissonModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/PoissonModel.R (rev 0)
+++ pkg/ROptEst/inst/scripts/PoissonModel.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,109 @@
+###############################################################################
+## Example: Poisson Family
+###############################################################################
+require(ROptEst)
+
+distroptions("TruncQuantile", 1e-10) # increases numerical support of Pois;
+ # i.e., increases precision of the
+ # computations
+## generates Poisson Family with theta = 4.5
+P <- PoisFamily(lambda = 4.5)
+P # show P
+plot(P) # plot of Pois(lambda = 4.5) and L_2 derivative
+checkL2deriv(P)
+
+# classical optimal IC
+IC0 <- optIC(model = P, risk = asCov())
+IC0 # show IC
+checkIC(IC0)
+Risks(IC0)
+plot(IC0) # plot IC
+
+# L_2 family + infinitesimal neighborhood
+RobP1 <- InfRobModel(center = P, neighbor = ContNeighborhood(radius = 0.5))
+RobP1 # show RobP1
+(RobP2 <- InfRobModel(center = P, neighbor = TotalVarNeighborhood(radius = 0.5)))
+
+## lower case radius
+lowerCaseRadius(L2Fam = P, ContNeighborhood(radius = 0.5), risk = asMSE())
+lowerCaseRadius(L2Fam = P, TotalVarNeighborhood(radius = 0.5), risk = asMSE())
+
+# MSE solution
+(IC1 <- optIC(model=RobP1, risk=asMSE()))
+checkIC(IC1)
+Risks(IC1)
+plot(IC1)
+
+(IC2 <- optIC(model=RobP2, risk=asMSE()))
+checkIC(IC2)
+Risks(IC2)
+plot(IC2)
+
+
+# lower case solutions
+(IC3 <- optIC(model=RobP1, risk=asBias()))
+checkIC(IC3)
+Risks(IC3)
+plot(IC3)
+
+(IC4 <- optIC(model=RobP2, risk=asBias()))
+checkIC(IC4)
+Risks(IC4)
+plot(IC4)
+
+# Hampel solution
+(IC5 <- optIC(model=RobP1, risk=asHampel(bound=clip(IC1))))
+checkIC(IC5)
+Risks(IC5)
+plot(IC5)
+
+(IC6 <- optIC(model=RobP2, risk=asHampel(bound=Risks(IC2)$asBias), maxiter = 200))
+checkIC(IC6)
+Risks(IC6)
+plot(IC6)
+
+
+# radius minimax IC
+(IC7 <- radiusMinimaxIC(L2Fam=P, neighbor=ContNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=0.5))
+checkIC(IC7)
+Risks(IC7)
+plot(IC7)
+
+(IC8 <- radiusMinimaxIC(L2Fam=P, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), loRad=0, upRad=Inf))
+checkIC(IC8)
+Risks(IC8)
+plot(IC8)
+
+# least favorable radius
+# (may take quite some time!)
+(r.rho1 <- leastFavorableRadius(L2Fam=P, neighbor=ContNeighborhood(),
+ risk=asMSE(), rho=0.5))
+(r.rho2 <- leastFavorableRadius(L2Fam=P, neighbor=TotalVarNeighborhood(),
+ risk=asMSE(), rho=1/3))
+
+## one-step estimation
+## Example: Rutherford-Geiger (1910)
+## cf. Feller~(1968), Section VI.7 (a)
+x <- c(rep(0, 57), rep(1, 203), rep(2, 383), rep(3, 525), rep(4, 532),
+ rep(5, 408), rep(6, 273), rep(7, 139), rep(8, 45), rep(9, 27),
+ rep(10, 10), rep(11, 4), rep(12, 0), rep(13, 1), rep(14, 1))
+
+## 0. mean (classical optimal)
+(est0 <- mean(x))
+
+## 1. Kolmogorov(-Smirnov) minimum distance estimator
+(est1 <- ksEstimator(x=x, Pois()))
+
+## 2. one-step estimation: radius interval
+## 2.1 small amount of contamination < 2%
+IC9 <- radiusMinimaxIC(L2Fam=PoisFamily(lambda=est1$lambda),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=1)
+(est21 <- oneStepEstimator(x, IC=IC9, start=est1$lambda))
+## 2.2 amount of contamination unknown
+IC10 <- radiusMinimaxIC(L2Fam=PoisFamily(lambda=est1$lambda),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(est22 <- oneStepEstimator(x, IC=IC10, start=est1$lambda))
+
+distroptions("TruncQuantile", 1e-5) # default
Added: pkg/ROptEst/inst/scripts/UnderOverShootRisk.R
===================================================================
--- pkg/ROptEst/inst/scripts/UnderOverShootRisk.R (rev 0)
+++ pkg/ROptEst/inst/scripts/UnderOverShootRisk.R 2008-02-19 06:09:45 UTC (rev 46)
@@ -0,0 +1,123 @@
+###############################################################################
+## Example: Normal Location
+###############################################################################
+system.time(require(ROptEst))
+
+## generates Normal Location Family with mean = 0
+N0 <- NormLocationFamily(mean=0)
+n <- 100
+tau <- qnorm(0.975)
+
+## classical optimal IC (radius = 0!)
+N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 0))
+N0.Rob2 <- InfRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = 0))
+
+system.time(IC0c <- optIC(model=N0.Rob1, risk=asUnOvShoot(width = tau)), gcFirst = TRUE)
+checkIC(IC0c)
+Risks(IC0c)
+system.time(IC0v <- optIC(model=N0.Rob2, risk=asUnOvShoot(width = tau)), gcFirst = TRUE)
+checkIC(IC0v)
+Risks(IC0v)
+
+## boundary case
+N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 2*tau*1/sqrt(2*pi)))
+N0.Rob2 <- InfRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = tau*1/sqrt(2*pi)))
+
+system.time(IC0c <- optIC(model=N0.Rob1, risk=asUnOvShoot(width = tau)), gcFirst = TRUE)
+checkIC(IC0c)
+Risks(IC0c)
+system.time(IC0v <- optIC(model=N0.Rob2, risk=asUnOvShoot(width = tau)), gcFirst = TRUE)
+checkIC(IC0v)
+Risks(IC0v)
+
+
+# L_2 family + infinitesimal resp. fixed neighborhood
+rad <- 0.5
+N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = rad))
+N0.Rob2 <- InfRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = rad/2))
+N0.Rob3 <- FixRobModel(center = N0, neighbor = ContNeighborhood(radius = rad/sqrt(n)))
+N0.Rob4 <- FixRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = rad/2/sqrt(n)))
+
+# asUnOvShoot solution
+N0.IC1 <- optIC(model = N0.Rob1, risk = asUnOvShoot(width = tau))
+checkIC(N0.IC1)
+Risks(N0.IC1)
+plot(N0.IC1)
+
+N0.IC2 <- optIC(model = N0.Rob2, risk = asUnOvShoot(width = tau))
+checkIC(N0.IC2)
+Risks(N0.IC2)
+plot(N0.IC2)
+
+# fiUnOvShoot solution
+N0.IC3 <- optIC(model=N0.Rob3, risk=fiUnOvShoot(width = tau/sqrt(n)), sampleSize = n)
+checkIC(N0.IC3)
+Risks(N0.IC3)
+plot(N0.IC3)
+
+N0.IC4 <- optIC(model=N0.Rob4, risk=fiUnOvShoot(width = tau/sqrt(n)), sampleSize = n)
+checkIC(N0.IC4)
+Risks(N0.IC4)
+plot(N0.IC4)
+
+# O(n^(-0.5))-corrected solution
+# in case of contamination neighborhoods
+N0.IC5 <- N0.IC1
+clipUp1 <- clipUp(N0.IC1)/as.vector(stand(N0.IC1))
+clipUp5 <- max(0, clipUp1 - rad*(rad + clipUp1*tau)/(sqrt(n)*2*tau*pnorm(-clipUp1)))
+stand5 <- 1/(2*pnorm(clipUp5)-1)
+clipUp(N0.IC5) <- stand5*clipUp5
+clipLo(N0.IC5) <- -clipUp(N0.IC5)
+stand(N0.IC5) <- as.matrix(stand5)
+Infos(N0.IC5) <- matrix(c("manually", "O(n^(-1/2))-corrected solution"), ncol=2,
+ dimnames=list(character(0), c("method", "message")))
+checkIC(N0.IC5)
+getRiskIC(N0.IC5, asUnOvShoot(width = tau), ContNeighborhood(radius=rad))
+getRiskIC(N0.IC5, fiUnOvShoot(width = tau/sqrt(n)), ContNeighborhood(radius=rad/sqrt(n)), sampleSize = n)
+
+# O(n^(-1))-corrected solution
+# in case of total variation neighborhoods
+N0.IC6 <- N0.IC2
+clipUp2 <- clipUp(N0.IC2)/as.vector(stand(N0.IC2))
+clipUp6 <- max(0, clipUp2 - tau*(2*clipUp2^2*rad/2 + tau*dnorm(clipUp2))/(6*n*pnorm(-clipUp2)))
+stand6 <- 1/(2*pnorm(clipUp6)-1)
+clipUp(N0.IC6) <- stand6*clipUp6
+clipLo(N0.IC6) <- -clipUp(N0.IC6)
+stand(N0.IC6) <- as.matrix(stand6)
+Infos(N0.IC6) <- matrix(c("manually", "O(n^(-1))-corrected solution"), ncol=2,
+ dimnames=list(character(0), c("method", "message")))
+checkIC(N0.IC6)
+getRiskIC(N0.IC6, asUnOvShoot(width = tau), TotalVarNeighborhood(radius=rad/2))
+getRiskIC(N0.IC6, fiUnOvShoot(width = tau/sqrt(n)), TotalVarNeighborhood(radius=rad/2/sqrt(n)), sampleSize = n)
+
+
+## estimation
+## 1. generate a contaminated sample
+ind <- rbinom(1e2, size = 1, prob = 0.05)
+X <- rnorm(1e2, mean = ind*4)
+summary(X)
+
+## 2. M estimation
+N0.Rob5 <- InfRobModel(center = NormLocationFamily(mean = 0),
+ neighbor = ContNeighborhood(radius = 0.5))
+N0.IC7 <- optIC(model=N0.Rob5, risk=asUnOvShoot(width = 1.960))
+(Mest1 <- locMEstimator(X, IC=N0.IC7))
+
+N0.Rob6 <- FixRobModel(center = NormLocationFamily(mean = 0),
+ neighbor = ContNeighborhood(radius = 0.05))
+N0.IC8 <- optIC(model = N0.Rob6, risk=fiUnOvShoot(width = 0.196), sampleSize = 1e2)
+(Mest2 <- locMEstimator(X, IC=N0.IC8))
+
+
+## 3. Kolmogorov(-Smirnov) minimum distance estimator
+(est0 <- ksEstimator(x=X, Norm(), param = "mean"))
+
+## 4. one-step estimation
+N0.Rob7 <- InfRobModel(center = NormLocationFamily(mean = est0$mean),
+ neighbor = ContNeighborhood(radius=0.5))
+N0.IC9 <- optIC(model=N0.Rob7, risk=asUnOvShoot(width = 1.960))
+(est1 <- oneStepEstimator(X, IC = N0.IC9, start = est0$mean))
+N0.Rob8 <- FixRobModel(center = NormLocationFamily(mean = est0$mean),
+ neighbor = ContNeighborhood(radius=0.05))
+N0.IC10 <- optIC(model=N0.Rob8, risk=fiUnOvShoot(width = 0.196), sampleSize = 1e2)
+(est2 <- oneStepEstimator(X, IC = N0.IC10, start = est0$mean))
More information about the Robast-commits
mailing list