Binomial Model
3.5 Implementation using R
θk(1−θ)m−k−0.5 l= 0,1, . . . , m−1 (3.4.1) If X ∼ Binom(m, θ) and Y ∼ Binom(m, θ0) with θ0 > θ, then Y is stochasti-cally strictly larger than X; i.e., P(Y > l)> P(X > l) for l = 0,1, . . . , m−1 . Thus, Fl(θ) is strictly decreasing in θ ∈ (0,1) for each l = 0,1, . . . , m−1 with
θ→0limFl(θ) = 0.5 and lim
θ→1Fl(θ) = −0.5 . Hence, Fl(θ) has a unique zero for each l= 0,1, . . . , m−1 . Furthermore, if θl∈(0,1) is the zero for some l, then 1−θl is the zero for m−1−l (by symmetry). In particular, we obtainθ0= 1−m√
0.5 (l= 0 ) and θm−1= m√
0.5 (l =m−1 ) are the roots of F0 and Fm−1 in (0,1) , respec-tively. In case m is odd, the median is also non-unique forθ= 0.5 (l= (m−1)/2 ).
But, in this case the distribution of Λ is symmetric to 0.5 and we can neglect the centering condition; confer Remark3.2.1(d).
In view of the previous considerations, the non-uniqueness of the centering con-stant zr is only a minor problem for the one-step construction of the corresponding optimally robust estimators and we may neglect this non-uniqueness in practice.
Since the minimum Kolmogorov(–Smirnov) distance estimator has the necessary properties (strict and √
n consistent on Uκ(θ) ⊃ Uv(θ) ⊃ Uc(θ) ) by Subsec-tion 6.3.2 ofRieder(1994) and is also well computable in this simple discrete model, we use it as our starting estimator.
3.5 Implementation using R
For a detailed description of theRpackageROptEst, which is part of ourRbundle RobASt, including the definition of the classes and methods mentioned below we refer to Appendix D.3. To generate a memberB of the binomial family with size m ∈ N and probability of success θ ∈(0,1) , we provide the generating function BinomFamily; i.e.,
> B <- BinomFamily(size = m, prob = θ)
One can also specify a transformation D ∈ R of the parameter θ by adding trafo=as.matrix(D)to the call of BinomFamily. The classical optimal (partial) IC (3.1.5) corresponding toBcan then be computed using
> IC0 <- optIC(model = B, risk = asCov())
102 Binomial Model
That is, calling the methodoptIC with an object of classL2ParamFamily and an object of classasCov returns an objectIC0of classIC.
For the computation of optimally robust ICs we provide the class InfRobModel which in addition to an object of class L2ParamFamilyincludes an unconditional infinitesimal neighborhood. The call
> RobB1 <- InfRobModel(center = B,
+ neighbor = ContNeighborhood(radius = r))
respectively
> RobB2 <- InfRobModel(center = B,
+ neighbor=TotalVarNeighborhood(radius = r))
generates an instanceRobB1, respectivelyRobB2of a binomial family with contam-ination, respectively total variation neighborhood and radius r ∈[0,∞] . In case we want to know, if the specified radius r is larger than the corresponding lower case radius ¯r, we can use
> lowerCaseRadius(L2Fam = B, neighbor = ContNeighborhood(),
+ risk = asMSE())
respectively
> lowerCaseRadius(L2Fam = B, neighbor = TotalVarNeighborhood(),
+ risk = asMSE())
The optimally robust IC with respect to the maximum asymptotic MSE (cf. (3.2.1), (3.2.17)) can then be computed via
> IC1 <- optIC(model = RobB1, risk = asMSE()) respectively
> IC2 <- optIC(model = RobB2, risk = asMSE())
That is, an object of class ContIC, respectively TotalVarIC is returned. In case r= 0 , the classical optimal IC and in case r=∞, the minimum bias solution (cf.
(3.2.5), (3.2.21)) are computed, respectively. The minimum bias solution can also be obtained by calling
> IC3 <- optIC(model = RobB1, risk = asBias()) respectively
> IC4 <- optIC(model = RobB2, risk = asBias())
In addition, it is also possible to determine the solutions to the corresponding Hampel type problem for a given clipping bound b∈ ωmin∗ , ω∗(ηh)
(∗=c, v) by using
> IC5 <- optIC(model = RobB1, risk = asHampel(bound = b))
3.5 Implementation usingR 103
respectively
> IC6 <- optIC(model = RobB2, risk = asHampel(bound = b))
In case b≥ω∗(ηh) , the classical optimal IC and in case b≤ω∗min, the minimum bias solution is returned, respectively.
If the radius is unknown, respectively unknown except to belong to some interval [a, b] with a∈[0,∞) and b∈(a,∞] , one can call
> IC7 <- radiusMinimaxIC(L2Fam = B, neighbor = ContNeighborhood(),
+ risk = asMSE(), loRad = a, upRad = b)
respectively
> IC8 <- radiusMinimaxIC(L2Fam = B,
+ neighbor = TotalVarNeighborhood(),
+ risk = asMSE(), loRad = a, upRad = b)
which computes the IC of the AL estimator which is radius–minimax in sense of Section 2.2. Finally, one can determine the least favorable radius rρ (ρ∈(0,1) ) and the corresponding MSE–inefficiency using
> r.rho1 <- leastFavorableRadius(L2Fam = B,
+ neighbor = ContNeighborhood(),
+ risk = asMSE(), rho = ρ)
respectively
> r.rho2 <- leastFavorableRadius(L2Fam = B,
+ neighbor = TotalVarNeighborhood(),
+ risk = asMSE(), rho = ρ)
which returns a list with membersrho,leastFavorableRadius and ineff. As pointed out in the following remark, the computation of the least favorable radius may take quite some time; in particular, in case of total variation neighborhoods.
Hence, intermediate results are printed to bridge this time gap.
Remark 3.5.1 The convergence of the algorithm for the computation of the opti-mally robust ICs is clearly faster in case of contamination neighborhoods. More pre-cisely, the computation time for the optimally robust IC is about 5−10 times larger in case of total variation neighborhoods. For instance, in case Binom (25,0.25) and
∗ = c, the computation time is about 0.25 seconds whereas in case ∗ = v it is about 1.65 seconds on an AMD Athlon with 2.5 GHz and 512 MB RAM usingR 2.0.1; conferR Development Core Team (2005). This, of course, also affects the computation of the radius–minimax IC (∗=c: approx. 2.1 sec. ∗ =v: approx.
33 sec.) and of the least favorable radius (∗=c: approx. 25 sec. ∗=v: approx.
8 min.). ////
104 Binomial Model
Having some (contaminated) sample X of Binom (m, θ) distributed data where m ∈ N is known, we can perform a robust estimation using the one-step con-struction. As starting estimator we provide the Kolmogorov(–Smirnov) minimum distance estimator which can be computed via
> est0 <- ksEstimator(x = X, distribution = Binom(size = m),
+ param = "prob")
With this initial estimator we then determine the corresponding optimally robust one-step estimator. If we are sure about the neighborhood radius r∈(0,∞) , we can use
> B0 <- BinomFamily(size = m, prob = est0)
> RobB3 <- InfRobModel(center = B0,
+ neighbor = ContNeighborhood(radius = r))
> IC9 <- optIC(model = RobB3, risk = asMSE())
> est1 <- oneStepEstimator(x = X, IC = IC9, start = est0)
respectively
> RobB4 <- InfRobModel(center = B0,
+ neighbor=TotalVarNeighborhood(radius = r))
> IC10 <- optIC(model = RobB4, risk = asMSE())
> est2 <- oneStepEstimator(x = X, IC = IC10, start = est0)
However, if the neighborhood radius is unknown except to lie in some interval [a, b]
(a∈[0,∞) , b∈(a,∞] ), we instead can proceed as follows
> IC11 <- radiusMinimaxIC(L2Fam = B0,
+ neighbor = ContNeighborhood(),
+ risk = asMSE(), loRad = a, upRad = b)
> est3 <- oneStepEstimator(x = X, IC = IC11, start = est0)
respectively
> IC12 <- radiusMinimaxIC(L2Fam = B0,
+ neighbor = TotalVarNeighborhood(),
+ risk = asMSE(), loRad = a, upRad = b)
> est4 <- oneStepEstimator(x = X, IC = IC12, start = est0)
Remark 3.5.2 After installing our R bundle RobASt one can find the R script BinomialModel.R, which contains some examples for the binomial model, in the directory “. . . /RHome/library/ROptEst/scripts/” where RHome stands for the
local home directory ofR. ////