Upper Record Values from a Univariate Normal Distribution
rec.normal.RdMaximum likelihood estimation of the two parameters of a univariate normal distribution when the observations are upper record values.
Usage
rec.normal(lmean = "identitylink", lsd = "loglink",
imean = NULL, isd = NULL, imethod = 1, zero = NULL)Arguments
- lmean, lsd
Link functions applied to the mean and sd parameters. See
Linksfor more choices.- imean, isd
Numeric. Optional initial values for the mean and sd. The default value
NULLmeans they are computed internally, with the help ofimethod.- imethod
Integer, either 1 or 2 or 3. Initial method, three algorithms are implemented. Choose the another value if convergence fails, or use
imeanand/orisd.- zero
Can be an integer vector, containing the value 1 or 2. If so, the mean or standard deviation respectively are modelled as an intercept only. Usually, setting
zero = 2will be used, if used at all. The default valueNULLmeans both linear/additive predictors are modelled as functions of the explanatory variables. SeeCommonVGAMffArgumentsfor more information.
Value
An object of class "vglmff"
(see vglmff-class).
The object is used by modelling functions
such as vglm,
and vgam.
References
Arnold, B. C. and Balakrishnan, N. and Nagaraja, H. N. (1998). Records, New York: John Wiley & Sons.
Note
This family function tries to solve a difficult problem,
and the larger the data set the better.
Convergence failure can commonly occur, and
convergence may be very slow,
so set maxit = 200, trace = TRUE, say.
Inputting good initial values are advised.
This family function uses the BFGS quasi-Newton update
formula for the
working weight matrices. Consequently the estimated
variance-covariance matrix may be inaccurate or
simply wrong! The
standard errors must be therefore treated with caution;
these are
computed in functions such as vcov()
and summary().
Examples
nn <- 10000; mymean <- 100
# First value is reference value or trivial record
Rdata <- data.frame(rawy = c(mymean, rnorm(nn, mymean, exp(3))))
# Keep only observations that are records:
rdata <- data.frame(y = unique(cummax(with(Rdata, rawy))))
fit <- vglm(y ~ 1, rec.normal, rdata, trace = TRUE, maxit = 200)
#> Iteration 1: loglikelihood = -30.664191
#> 12 weight matrices not updated out of 13
#> Iteration 2: loglikelihood = -26.76061
#> 11 weight matrices not updated out of 13
#> Iteration 3: loglikelihood = -26.743649
#> 12 weight matrices not updated out of 13
#> Iteration 4: loglikelihood = -26.738771
#> 12 weight matrices not updated out of 13
#> Iteration 5: loglikelihood = -26.738448
#> 12 weight matrices not updated out of 13
#> Iteration 6: loglikelihood = -26.738398
#> 12 weight matrices not updated out of 13
#> Iteration 7: loglikelihood = -26.738365
#> 12 weight matrices not updated out of 13
#> Iteration 8: loglikelihood = -26.738332
#> 12 weight matrices not updated out of 13
#> Iteration 9: loglikelihood = -26.7383
#> 12 weight matrices not updated out of 13
#> Iteration 10: loglikelihood = -26.738267
#> Iteration 11: loglikelihood = -26.635721
#> Iteration 12: loglikelihood = -26.633118
#> 1 weight matrices not updated out of 13
#> Iteration 13: loglikelihood = -26.632205
#> 1 weight matrices not updated out of 13
#> Iteration 14: loglikelihood = -26.631345
#> Iteration 15: loglikelihood = -26.605746
#> Iteration 16: loglikelihood = -26.605587
#> 3 weight matrices not updated out of 13
#> Iteration 17: loglikelihood = -26.605436
#> 12 weight matrices not updated out of 13
#> Iteration 18: loglikelihood = -26.605367
#> 1 weight matrices not updated out of 13
#> Iteration 19: loglikelihood = -26.605291
#> 1 weight matrices not updated out of 13
#> Iteration 20: loglikelihood = -26.605279
#> 1 weight matrices not updated out of 13
#> Iteration 21: loglikelihood = -26.605277
#> 1 weight matrices not updated out of 13
#> Iteration 22: loglikelihood = -26.605276
coef(fit, matrix = TRUE)
#> mean loglink(sd)
#> (Intercept) 92.58273 2.888791
Coef(fit)
#> mu sd
#> 92.58273 17.97157
summary(fit)
#>
#> Call:
#> vglm(formula = y ~ 1, family = rec.normal, data = rdata, trace = TRUE,
#> maxit = 200)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 92.5827 11.9673 7.736 1.02e-14 ***
#> (Intercept):2 2.8888 0.1794 16.101 < 2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: mean, loglink(sd)
#>
#> Log-likelihood: -26.6053 on 24 degrees of freedom
#>
#> Number of Fisher scoring iterations: 22
#>