This function calculates the log-likelihood of an independent and identically distributed (iid) sample from a distribution. See Details.
Value
If only the distr argument is supplied, ll() returns a function.
If both distr and x are supplied, ll() returns a numeric, the value of
the log-likelihood function.
Details
The log-likelihood functions are provided in two forms: the ll<name>
distribution-specific version that follows the stats package conventions, and
the S4 generic ll. Examples for the ll<name> version are included in the
distribution-specific help pages, e.g. ?Beta (all distributions can be
found in the See Also section of the Distributions help page).
As with the d(), p(), q(), r() methods, ll() can be supplied only
with distr to return the log-likelihood function (i.e. it can be used as a
functional), or with both distr and x to be evaluated directly.
In some distribution families like beta and gamma, the MLE cannot be explicitly derived and numerical optimization algorithms have to be employed. Even in “good" scenarios, with plenty of observations and a smooth optimization function, extra care should be taken to ensure a fast and right convergence if possible. Two important steps are taken in package in this direction:
The log-likelihood function is analytically calculated for each distribution family, so that constant terms with respect to the parameters can be removed, leaving only the sufficient statistics as a requirement for the function evaluation.
Multidimensional problems are reduced to unidimensional ones by utilizing the score equations.
The resulting function that is inserted in the optimization algorithm is
called lloptim(), and is not to be confused with the actual log-likelihood
function ll(). The corresponding derivative is called dlloptim(). These
functions are used internally and are not exported.
Therefore, whenever numerical computation of the MLE is required, optim()
is called to optimize lloptim(), using the ME or SAME as the starting point
(user's choice), and the L-BFGS-U optimization algorithm, with lower and
upper limits defined by default as the parameter space boundary. Illustrative
examples can be found in the package vignette.
Examples
# -----------------------------------------------------
# Beta Distribution Example
# -----------------------------------------------------
# Create the distribution
a <- 3
b <- 5
D <- Beta(a, b)
# ------------------
# dpqr Functions
# ------------------
d(D, c(0.3, 0.8, 0.5)) # density function
#> [1] 2.268945 0.107520 1.640625
p(D, c(0.3, 0.8, 0.5)) # distribution function
#> [1] 0.3529305 0.9953280 0.7734375
qn(D, c(0.4, 0.8)) # inverse distribution function
#> [1] 0.3205858 0.5167578
x <- r(D, 100) # random generator function
# alternative way to use the function
df <- d(D) ; df(x) # df is a function itself
#> [1] 1.7851371 1.6158350 1.8831781 2.3045040 2.2205275 1.5793396 1.1770878
#> [8] 2.3042618 1.2800954 0.6658872 2.0576716 2.2848265 1.6149912 1.0887630
#> [15] 2.2317215 1.6721279 2.1187735 0.2271411 1.7116001 1.8887920 2.1966674
#> [22] 0.6208597 0.4909717 1.9340585 2.2601533 0.9333753 2.2851968 1.6134685
#> [29] 1.8938827 2.2891671 2.3039898 1.6144085 2.2729664 1.5983859 1.6480103
#> [36] 1.9552379 1.3160884 0.4814330 1.7809212 1.3228666 1.3184469 2.0666304
#> [43] 2.2827307 2.0837601 1.5420218 2.1792920 1.2426316 2.1115211 1.9887922
#> [50] 1.7716759 2.2223997 1.7971389 1.4612738 0.5583837 2.1420633 2.0977546
#> [57] 1.3226088 1.3105294 0.8628855 1.7648335 1.7796204 2.0334251 1.5808053
#> [64] 1.4721859 1.9030829 1.0340463 1.4202481 1.7353987 2.2324880 1.4847137
#> [71] 2.2580652 1.1813369 1.6089141 1.5267540 0.6130996 2.2864078 2.1086791
#> [78] 2.2795250 1.4493381 2.1531247 1.0302573 2.2977963 2.3042097 2.3023549
#> [85] 2.2933047 0.1138212 2.2778024 1.9848331 2.2984321 2.2476951 0.7530797
#> [92] 2.2385515 1.1660503 2.0617423 1.3656245 1.9102785 2.1218297 2.2878026
#> [99] 1.5673847 1.6756807
# ------------------
# Moments
# ------------------
mean(D) # Expectation
#> [1] 0.375
var(D) # Variance
#> [1] 0.02604167
sd(D) # Standard Deviation
#> [1] 0.1613743
skew(D) # Skewness
#> [1] 0.3098387
kurt(D) # Excess Kurtosis
#> [1] 0.04
entro(D) # Entropy
#> [1] -0.4301508
finf(D) # Fisher Information Matrix
#> shape1 shape2
#> shape1 0.2617971 -0.13313701
#> shape2 -0.1331370 0.08818594
# List of all available moments
mom <- moments(D)
mom$mean # expectation
#> [1] 0.375
# ------------------
# Point Estimation
# ------------------
ll(D, x)
#> [1] 45.17075
llbeta(x, a, b)
#> [1] 45.17075
ebeta(x, type = "mle")
#> $shape1
#> [1] 3.154097
#>
#> $shape2
#> [1] 5.279473
#>
ebeta(x, type = "me")
#> $shape1
#> [1] 3.015129
#>
#> $shape2
#> [1] 5.065686
#>
ebeta(x, type = "same")
#> $shape1
#> [1] 3.107776
#>
#> $shape2
#> [1] 5.221341
#>
mle(D, x)
#> $shape1
#> [1] 3.154097
#>
#> $shape2
#> [1] 5.279473
#>
me(D, x)
#> $shape1
#> [1] 3.015129
#>
#> $shape2
#> [1] 5.065686
#>
same(D, x)
#> $shape1
#> [1] 3.107776
#>
#> $shape2
#> [1] 5.221341
#>
e(D, x, type = "mle")
#> $shape1
#> [1] 3.154097
#>
#> $shape2
#> [1] 5.279473
#>
mle("beta", x) # the distr argument can be a character
#> $shape1
#> [1] 3.154097
#>
#> $shape2
#> [1] 5.279473
#>
# ------------------
# Estimator Variance
# ------------------
vbeta(a, b, type = "mle")
#> shape1 shape2
#> shape1 16.44844 24.83272
#> shape2 24.83272 48.83039
vbeta(a, b, type = "me")
#> shape1 shape2
#> shape1 17.64848 26.56970
#> shape2 26.56970 51.39394
vbeta(a, b, type = "same")
#> shape1 shape2
#> shape1 16.57719 24.96198
#> shape2 24.96198 49.01071
avar_mle(D)
#> shape1 shape2
#> shape1 16.44844 24.83272
#> shape2 24.83272 48.83039
avar_me(D)
#> shape1 shape2
#> shape1 17.64848 26.56970
#> shape2 26.56970 51.39394
avar_same(D)
#> shape1 shape2
#> shape1 16.57719 24.96198
#> shape2 24.96198 49.01071
v(D, type = "mle")
#> shape1 shape2
#> shape1 16.44844 24.83272
#> shape2 24.83272 48.83039