# Repeat the simulations described in Harter and Moore (1966)
SIMS <- 100 # Number of simulations (change this to 1000)
mu.save <- sd.save <- rep(NA, len = SIMS)
r1 <- 0; r2 <- 4; nn <- 20
for(sim in 1:SIMS) {
y <- sort(rnorm(nn))
y <- y[(1+r1):(nn-r2)] # Delete r1 smallest and r2 largest
fit <- vglm(y ~ 1, dcennormal1(r1 = r1, r2 = r2))
mu.save[sim] <- predict(fit)[1,1]
sd.save[sim] <- exp(predict(fit)[1,2]) # Assumes a log link and ~ 1
}
c(mean(mu.save), mean(sd.save)) # Should be c(0,1)
c(sd(mu.save), sd(sd.save))
# Data from Sarhan and Greenberg (1962); MLEs are mu = 9.2606, sd = 1.3754
strontium90 <- data.frame(y = c(8.2, 8.4, 9.1, 9.8, 9.9))
fit <- vglm(y ~ 1, dcennormal1(r1 = 2, r2 = 3, isd = 6), strontium90, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
Run the code above in your browser using DataLab