leipnik {VGAM} | R Documentation |
Estimates the two parameters of a (transformed) Leipnik distribution by maximum likelihood estimation.
leipnik(lmu = "logit", llambda = "loge", emu=list(), elambda=list(), imu = NULL, ilambda = NULL)
lmu, llambda |
Link function for the mu and lambda parameters.
See Links for more choices.
|
imu, ilambda |
Numeric. Optional initial values for mu and
lambda.
|
emu, elambda |
List. Extra argument for each of the links.
See earg in Links for general information.
|
The (transformed) Leipnik distribution has density function
f(y;mu,lambda) = (y(1-y))^(-1/2) * (1 + (y-mu)^2 / (y*(1-y)))^(-lambda/2) / Beta((lambda+1)/2, 1/2)
where 0 < y < 1 and lambda > -1. The mean is mu (returned as the fitted values) and the variance is 1/lambda.
Jorgensen (1997) calls the above the transformed Leipnik distribution, and if y = (x+1)/2 and mu = (theta+1)/2, then the distribution of X as a function of x and theta is known as the the (untransformed) Leipnik distribution. Here, both x and theta are in (-1,1).
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions such as vglm
,
rrvglm
and vgam
.
If llambda="identity"
then it is possible that the
lambda
estimate becomes less than -1, i.e., out of
bounds. One way to stop this is to choose llambda="loge"
,
however, lambda
is then constrained to be positive.
Convergence may be slow or fail.
Until better initial value estimates are forthcoming try assigning the
argument ilambda
some numerical value if it fails to converge.
Currently, Newton-Raphson is implemented, not Fisher scoring.
Currently, this family function probably only really works for
intercept-only models, i.e., y ~ 1
in the formula.
T. W. Yee
Jorgensen, B. (1997) The Theory of Dispersion Models. London: Chapman & Hall
Johnson, N. L. and Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions, 2nd edition, Volume 2, New York: Wiley. (pages 612–617).
y = rnorm(n=2000, mean=0.5, sd=0.1) # Not good data fit = vglm(y ~ 1, leipnik(ilambda=1), tr=TRUE, checkwz=FALSE) fit = vglm(y ~ 1, leipnik(ilambda=1,llam=logoff, elam=list(offset=1)), trace=TRUE, cri="coef") fitted(fit)[1:5] mean(y) summary(fit) coef(fit, matrix=TRUE) Coef(fit) sum(weights(fit)) # sum of the prior weights sum(weights(fit, type="w")) # sum of the working weights