Zeta Distribution Family Function

Estimates the parameter of the zeta distribution.

Usage

zetaff(lshape = "loglink", ishape = NULL, gshape = 1 + exp(-seq(7)),
       zero = NULL)

Arguments

lshape, ishape, zero

These arguments apply to the (positive) parameter $p$. See Links for more choices. Choosing loglog constrains $p>1$, but may fail if the maximum likelihood estimate is less than one. See CommonVGAMffArguments for more information.

gshape

See CommonVGAMffArguments for more information.

Details

In this long tailed distribution the response must be a positive integer. The probability function for a response $Y$ is $$P(Y=y) = 1/[y^{p+1} \zeta(p+1)],\ \ \ p>0,\ \ \ y=1,2,...$$ where $\zeta$ is Riemann's zeta function. The parameter $p$ is positive, therefore a log link is the default. The mean of $Y$ is $\mu = \zeta(p) / \zeta(p+1)$ (provided $p>1$) and these are the fitted values. The variance of $Y$ is $\zeta(p-1) / \zeta(p+1) - \mu^2$ provided $p>2$.

It appears that good initial values are needed for successful convergence. If convergence is not obtained, try several values ranging from values near 0 to values about 10 or more.

Multiple responses are handled.

Value

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

References

pp.527– of Chapter 11 of Johnson N. L., Kemp, A. W. and Kotz S. (2005). Univariate Discrete Distributions, 3rd edition, Hoboken, New Jersey: Wiley.

Knight, K. (2000). Mathematical Statistics. Boca Raton, FL, USA: Chapman & Hall/CRC Press.

Author

T. W. Yee

Note

The zeta function may be used to compute values of the zeta function.

Examples

zdata <- data.frame(y = 1:5, w =  c(63, 14, 5, 1, 2))  # Knight, p.304
fit <- vglm(y ~ 1, zetaff, data = zdata, trace = TRUE, weight = w, crit = "c")
#> Iteration 1: coefficients = 0.49363036
#> Iteration 2: coefficients = 0.51981018
#> Iteration 3: coefficients = 0.5203144
#> Iteration 4: coefficients = 0.52031458
#> Iteration 5: coefficients = 0.52031458
(phat <- Coef(fit))  # 1.682557
#>    shape 
#> 1.682557 
with(zdata, cbind(round(dzeta(y, phat) * sum(w), 1), w))
#>            w
#> [1,] 66.4 63
#> [2,] 10.3 14
#> [3,]  3.5  5
#> [4,]  1.6  1
#> [5,]  0.9  2

with(zdata, weighted.mean(y, w))
#> [1] 1.411765
fitted(fit, matrix = FALSE)
#>       [,1]
#> 1 1.633302
#> 2 1.633302
#> 3 1.633302
#> 4 1.633302
#> 5 1.633302
predict(fit)
#>      loglink(shape)
#> [1,]      0.5203146
#> [2,]      0.5203146
#> [3,]      0.5203146
#> [4,]      0.5203146
#> [5,]      0.5203146

# The following should be zero at the MLE:
with(zdata, mean(log(rep(y, w))) + zeta(1+phat, deriv = 1) / zeta(1+phat))
#>         shape 
#> -2.775558e-17