1
2# This library is free software; you can redistribute it and/or
3# modify it under the terms of the GNU Library General Public
4# License as published by the Free Software Foundation; either
5# version 2 of the License, or (at your option) any later version.
6#
7# This library is distributed in the hope that it will be useful,
8# but WITHOUT ANY WARRANTY; without even the implied warranty of
9# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10# GNU Library General Public License for more details.
11#
12# You should have received a copy of the GNU Library General
13# Public License along with this library; if not, write to the
14# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
15# MA  02111-1307  USA
16
17
18################################################################################
19# FUNCTION:            GENERALIZED DISTRIBUTION:
20#  hypFit               Fits parameters of a hyperbolic density
21################################################################################
22
23
24hypFit <-
25function(x, alpha = 1, beta = 0, delta = 1, mu = 0, scale = TRUE,
26    doplot = TRUE, span = "auto", trace = TRUE,
27    title = NULL, description = NULL, ...)
28{
29    # A function implemented by Diethelm Wuertz
30
31    # Description:
32    #   Fits parameters of a hyperbolic density
33
34    # FUNCTION:
35
36    # Transform:
37    x.orig = x
38    x = as.vector(x)
39    if (scale) {
40        SD = sd(x)
41        x = x / SD
42    }
43
44    # Settings:
45    CALL = match.call()
46
47    # Log-likelihood Function:
48    ehypmle = function(x, y = x, trace) {
49        if (abs(x[2]) >= x[1]) return(1e99)
50        f = -sum(log(dhyp(y, x[1], x[2], x[3], x[4])))
51        # Print Iteration Path:
52        if (trace) {
53            cat("\n Objective Function Value:  ", -f)
54            cat("\n Parameter Estimates:       ", x[1], x[2], x[3], x[4], "\n")
55        }
56        f
57    }
58
59    # Minimization:
60    eps = 1e-10
61    BIG = 1000
62    r = nlminb(start = c(alpha, beta, delta, mu), objective = ehypmle,
63        lower = c(eps, -BIG, eps, -BIG), upper = BIG, y = x, trace = trace)
64    names(r$par) <- c("alpha", "beta", "delta", "mu")
65
66    # Add Title and Description:
67    if (is.null(title)) title = "Hyperbolic Parameter Estimation"
68    if (is.null(description)) description = description()
69
70    # Result:
71    if (scale) {
72        r$par = r$par / c(SD, SD, 1/SD, 1/SD)
73        r$objective = ehypmle(r$par, y = as.vector(x.orig), trace = trace)
74    }
75    fit = list(estimate = r$par, minimum = -r$objective, code = r$convergence)
76
77    # Optional Plot:
78    if (doplot) {
79        x = as.vector(x.orig)
80        if (span == "auto") span = seq(min(x), max(x), length = 51)
81        z = density(x, n = 100, ...)
82        x = z$x[z$y > 0]
83        y = z$y[z$y > 0]
84        y.points = dhyp(span, r$par[1], r$par[2], r$par[3], r$par[4])
85        ylim = log(c(min(y.points), max(y.points)))
86        plot(x, log(y), xlim = c(span[1], span[length(span)]),
87            ylim = ylim, type = "p", xlab = "x", ylab = "log f(x)", ...)
88        title("HYP Parameter Estimation")
89        lines(x = span, y = log(y.points), col = "steelblue")
90    }
91
92
93    # Return Value:
94    new("fDISTFIT",
95        call = as.call(CALL),
96        model = "Hyperbolic Distribution",
97        data = as.data.frame(x.orig),
98        fit = fit,
99        title = title,
100        description = description )
101}
102
103
104################################################################################
105
106