Univariate Normal Distribution with Double Censoring

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

double.cens.normal(r1 = 0, r2 = 0, lmu = "identitylink", lsd =
       "loglink", imu = NULL, isd = NULL, zero = "sd")

Arguments

r1, r2

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

lmu, lsd

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

imu, isd, zero

See CommonVGAMffArguments for more information.

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.

Value

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

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.

Author

T. W. Yee

Note

This family function only handles a vector or one-column matrix response. The weights argument, if used, are interpreted as frequencies, therefore it must be a vector with positive integer values.

With no censoring at all (the default), it is better (and equivalent) to use uninormal.

See also

uninormal, cens.normal, tobit.

Examples

if (FALSE)  # Repeat the simulations of Harter & Moore (1966)
SIMS <- 100  # Number of simulations (change this to 1000)
mu.save <- sd.save <- rep(NA, len = SIMS)
#> Error: object 'SIMS' not found
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, double.cens.normal(r1 = r1, r2 = r2))
  mu.save[sim] <- predict(fit)[1, 1]
  sd.save[sim] <- exp(predict(fit)[1, 2])  # Assumes a log link & ~ 1
}
#> Error: object 'SIMS' not found
c(mean(mu.save), mean(sd.save))  # Should be c(0,1)
#> Error: object 'mu.save' not found
c(sd(mu.save), sd(sd.save))
#> Error: object 'mu.save' not found
 # \dontrun{}

# Data from Sarhan & 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, double.cens.normal(r1 = 2, r2 = 3, isd = 6),
            data = strontium90, trace = TRUE)
#> Iteration 1: loglikelihood = -15.69275
#> Iteration 2: loglikelihood = -14.25529
#> Iteration 3: loglikelihood = -13.50897
#> Iteration 4: loglikelihood = -13.32076
#> Iteration 5: loglikelihood = -13.30763
#> Iteration 6: loglikelihood = -13.3075
#> Iteration 7: loglikelihood = -13.3075
#> Iteration 8: loglikelihood = -13.3075
coef(fit, matrix = TRUE)
#>                   mu loglink(sd)
#> (Intercept) 9.260564   0.3187177
Coef(fit)
#>       mu       sd 
#> 9.260564 1.375363