Statistics 5601 (Geyer, Spring 2006) Examples: Sandwich Estimator

General Instructions
Theory
The Cauchy Location Problem
The Cauchy Location Problem with Sandwich
The Cauchy Location-Scale Problem with Sandwich

General Instructions

To do each example, just click the "Submit" button. You do not have to type in any R instructions or specify a dataset. That's already done for you.

Theory

This page goes with notes, the PDF form of which is found on the theory page.

We just repeat the R computations in that note in Rweb form here.

The Cauchy Location Problem with Sandwich

logl <- expression(- log(1 + (x - theta)^2)) loglfun <- deriv3(logl, "theta", c("theta", "x")) loglfun theta.start <- median(x) ##### does work vectorwise ##### loglfun(theta.start, x = x) fred <- function(theta, x) { sally <- loglfun(theta, x) result <- sum(- sally) attr(result, "gradient") <- sum(- attr(sally, "gradient")) return(result) } out <- nlm(fred, median(x), x = x) print(out$estimate) print(out$gradient) vfun <- function(theta, x) { sally <- loglfun(theta, x) sum(attr(sally, "gradient")^2) } jfun <- function(theta, x) { sally <- loglfun(theta, x) sum(- attr(sally, "hessian")) } theta.hat <- out$estimate Vhat <- vfun(theta.hat, x) Jhat <- jfun(theta.hat, x) print(Vhat) print(Jhat) conf.level <- 0.95 crit <- qnorm((1 + conf.level) / 2) theta.hat + crit * c(-1, 1) * sqrt(Vhat / Jhat^2)

External Data Entry

Enter a dataset URL :

For security reasons, Rweb forbids the use of the eval function. This makes the code in the form above much messier than the code in the notes, but both calculate the same things.

The Cauchy Location-Scale Problem with Sandwich

logl <- expression(- log(1 + ((x - mu) / sigma)^2) - log(sigma)) loglfun <- deriv3(logl, c("mu", "sigma"), c("mu", "sigma", "x")) loglfun mu.start <- median(x) sigma.start <- IQR(x) / 2 ##### does work vectorwise ##### loglfun(mu.start, sigma.start, x = x) fred <- function(theta, x) { mu <- theta[1] sigma <- theta[2] sally <- loglfun(mu, sigma, x) result <- sum(- sally) attr(result, "gradient") <- apply(- attr(sally, "gradient"), 2, sum) return(result) } theta.start <- c(mu.start, sigma.start) out <- nlm(fred, theta.start, x = x) print(out$estimate) print(out$gradient) vfun <- function(theta, x) { mu <- theta[1] sigma <- theta[2] sally <- loglfun(mu, sigma, x) herman <- attr(sally, "gradient") t(herman) %*% herman } jfun <- function(theta, x) { mu <- theta[1] sigma <- theta[2] sally <- loglfun(mu, sigma, x) apply(- attr(sally, "hessian"), c(2, 3), sum) } theta.hat <- out$estimate Vhat <- vfun(theta.hat, x) Jhat <- jfun(theta.hat, x) Mhat <- solve(Jhat) %*% Vhat %*% solve(Jhat) print(Vhat) print(Jhat) print(Mhat) conf.level <- 0.95 crit <- qnorm((1 + conf.level) / 2) theta.hat[1] + crit * c(-1, 1) * sqrt(Mhat[1, 1]) theta.hat[2] + crit * c(-1, 1) * sqrt(Mhat[2, 2])