inst/scripts/LognormalAndNormalModel.R

###############################################################################
## Example: Lognormal Scale and Normal Location
###############################################################################
require(ROptEst)
options("newDevice"=TRUE)

## generates Lognormal Scale Family with meanlog = 0, sdlog = 1
LN1 <- LnormScaleFamily() 
LN1        # show LN1
plot(LN1)  # plot of Lnorm() 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$value))
checkIC(LN1.IC6)
Risks(LN1.IC6)
plot(LN1.IC6)

N0.IC6 <- optIC(model=N0.Rob2, risk=asHampel(bound=Risks(N0.IC2)$asBias$value))
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, loRad0=2e-3))
checkIC(LN1.IC7)
Risks(LN1.IC7)
plot(LN1.IC7)

(N0.IC7 <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(), 
                risk=asMSE(), loRad=0, upRad=0.5, loRad0=2e-3))
checkIC(N0.IC7)
Risks(N0.IC7)
plot(N0.IC7)

(LN1.IC8 <- radiusMinimaxIC(L2Fam=LN1, neighbor=TotalVarNeighborhood(), 
                risk=asMSE(), loRad=0, upRad=0.25, loRad0=0.04))
checkIC(LN1.IC8)
Risks(LN1.IC8)
plot(LN1.IC8)

(N0.IC8 <- radiusMinimaxIC(L2Fam=N0, neighbor=TotalVarNeighborhood(), 
                risk=asMSE(), loRad=0, upRad=0.25, loRad0=0.04))
checkIC(N0.IC8)
Risks(N0.IC8)
plot(N0.IC8)


## least favorable radius
(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))

## For estimation use function roptest
ind <- rbinom(1e2, size=1, prob=0.05) 
x <- rnorm(1e2, mean=(1-ind)+ind*9)
y <- exp(x)

## 1-step: contamination known
est1 <- roptest(x, eps = 0.05, L2Fam = NormLocationFamily())
est2 <- roptest(y, eps = 0.05, L2Fam = LnormScaleFamily())

## k-step: contamination known
est3 <- roptest(x, eps = 0.05, L2Fam = NormLocationFamily(), steps = 3)
est4 <- roptest(y, eps = 0.05, L2Fam = LnormScaleFamily(), steps = 3)

## comparison
estimate(est1)
log(estimate(est2))
estimate(est3)
log(estimate(est4))

## confidence intervals
confint(est1, symmetricBias())
confint(est2, symmetricBias())
confint(est3, symmetricBias())
confint(est4, symmetricBias())

Try the ROptEst package in your browser

Any scripts or data that you put into this service are public.

ROptEst documentation built on April 6, 2019, 3:01 a.m.