A test of Stein estimator

Stein.Rmd

---
title: "Stein estimator"
author: "ITÔ, Hiroki"
date: "2018/3/22"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Define sum of squared error

```{r sse}
sse <- function(theta.star, theta) {
  sum((theta.star - theta)^2)
}
```

## Define Stein estimator
```{r stein}
stein <- function(y, sigma) {
  n <- length(y)
  if (n <= 3)
    stop("n must be larger than 3.")
  s2 <- sum((y - mean(y))^2) / (n - 3)
  a <- sigma^2 / s2
  theta.stein <- (1 - a) * y + a * mean(y)
}
```

## Simulation

### Sample from uniform distribution

```{r sim1}
set.seed(123)
n <- 64
theta <- runif(n, 0, 128)
sigma <- 2
sse1 <- sapply(seq_len(2000), function(i) {
  y <- rnorm(n, theta, sigma)
  theta.stein <- stein(y, sigma)
  c(sse(y, theta), sse(theta.stein, theta))
})
apply(sse1, 1, mean)
```

### Sample from normal distribution

```{r sim2}
set.seed(123)
n <- 64
theta <- rnorm(n, 0, 1)
sigma <- 2
sse2 <- sapply(seq_len(2000), function(i) {
  y <- rnorm(n, theta, sigma)
  theta.stein <- stein(y, sigma)
  c(sse(y, theta), sse(theta.stein, theta))
})
apply(sse2, 1, mean)
```

### Sample from uniform distribution with a long range

```{r sim3}
set.seed(123)
n <- 64
theta <- runif(n, -1e+5, 1e+5)
sigma <- 1
sse3 <- sapply(1:2000, function(i) {
  y <- rnorm(n, theta, sigma)
  theta.stein <- stein(y, sigma)
  c(sse(y, theta), sse(theta.stein, theta))
})
apply(sse3, 1, mean)
```