dcennormal1: Univariate Normal Distribution with Double Censoring

Description

Maximum likelihood estimation of the two parameters of a univariate normal distribution when there is double censoring.

Usage

dcennormal1(r1 = 0, r2 = 0, lmu = "identity", lsd = "loge",
            emu = list(), esd = list(),
            imu = NULL, isd = NULL, zero = 2)

Arguments

r1, r2

Integers. Number of smallest and largest values censored, respectively.

lmu, lsd, emu, esd

Parameter link functions and its extra arguments applied to the mean and standard deviation. See Links for more choices.

imu, isd, zero

See CommonVGAMffArguments for more information.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Details

This family function uses the Fisher information matrix given in Harter and Moore (1966). The matrix is not diagonal if either r1 or r2 are positive.

By default, the mean is the first linear/additive predictor and the log of the standard deviation is the second linear/additive predictor.

References

Harter, H. L. and Moore, A. H. (1966) Iterative maximum-likelihood estimation of the parameters of normal populations from singly and doubly censored samples. Biometrika, 53, 205--213.

Examples

Run this code

# 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

Last chance! 50% off unlimited learning