Calculate a generalized logarithmic mean / extended mean.
Value
generalized_logmean() and extended_mean() return a
function:
function(a, b, tol = .Machine$double.eps^0.5){...}This computes the component-wise generalized logarithmic mean of order
r, or the extended mean of orders r and s, of a and b.
logmean() returns a numeric vector, the same length as
max(length(a), length(b)), giving the component-wise logarithmic mean
of a and b.
Details
The function extended_mean() returns a function to compute the
component-wise extended mean of a and b of orders r and
s. See Bullen (2003, p. 393) for a definition. This is also called
the difference mean, Stolarsky mean, or extended mean-value mean.
The function generalized_logmean() returns a function to compute the
component-wise generalized logarithmic mean of a and b of
order r. See Bullen (2003, p. 385) for a definition, or
https://en.wikipedia.org/wiki/Stolarsky_mean. The generalized
logarithmic mean is a special case of the extended mean, corresponding to
extended_mean(r, 1)(), but is more commonly used for price indexes.
The function logmean() returns the ordinary component-wise
logarithmic mean of a and b, and corresponds to
generalized_logmean(1)().
Both a and b should be strictly positive. This is not
enforced, but the results may not make sense when the generalized
logarithmic mean / extended mean is not defined. The usual recycling rules
apply when a and b are not the same length.
By definition, the generalized logarithmic mean / extended mean of a
and b is a when a == b. The tol argument is used
to test equality by checking if abs(a - b) <= tol. The default value
is the same as all.equal(). Setting tol = 0
tests for exact equality, but can give misleading results when a and
b are computed values. In some cases it's useful to multiply
tol by a scale factor, such as max(abs(a), abs(b)). This often
doesn't matter when making price indexes, however, as a and b
are usually around 1.
Note
generalized_logmean() can be defined on the extended real line,
so that r = -Inf / Inf returns pmin()/pmax(), to agree with the
definition in, e.g., Bullen (2003). This is not implemented, and r
must be finite as in the original formulation by Stolarsky (1975).
References
Balk, B. M. (2008). Price and Quantity Index Numbers. Cambridge University Press.
Bullen, P. S. (2003). Handbook of Means and Their Inequalities. Springer Science+Business Media.
Stolarsky, K. B. (1975). Generalizations of the Logarithmic Mean. Mathematics Magazine, 48(2): 87-92.
See also
transmute_weights() uses the extended mean to turn a generalized
mean of order \(r\) into a generalized mean of order \(s\).
Other means:
generalized_mean(),
lehmer_mean(),
nested_mean()
Examples
x <- 8:5
y <- 1:4
#---- Comparing logarithmic means and generalized means ----
# The arithmetic and geometric means are special cases of the
# generalized logarithmic mean
all.equal(generalized_logmean(2)(x, y), (x + y) / 2)
#> [1] TRUE
all.equal(generalized_logmean(-1)(x, y), sqrt(x * y))
#> [1] TRUE
# The logarithmic mean lies between the arithmetic and geometric means
# because the generalized logarithmic mean is increasing in r
all(logmean(x, y) < (x + y) / 2) &
all(logmean(x, y) > sqrt(x * y))
#> [1] TRUE
# The harmonic mean cannot be expressed as a logarithmic mean, but can
# be expressed as an extended mean
all.equal(extended_mean(-2, -1)(x, y), 2 / (1 / x + 1 / y))
#> [1] TRUE
# The quadratic mean is also a type of extended mean
all.equal(extended_mean(2, 4)(x, y), sqrt(x^2 / 2 + y^2 / 2))
#> [1] TRUE
# As are heronian and centroidal means
all.equal(
extended_mean(0.5, 1.5)(x, y),
(x + sqrt(x * y) + y) / 3
)
#> [1] TRUE
all.equal(
extended_mean(2, 3)(x, y),
2 / 3 * (x^2 + x * y + y^2) / (x + y)
)
#> [1] TRUE
#---- Approximating the logarithmic mean ----
# The logarithmic mean can be approximated as a convex combination of
# the arithmetic and geometric means that gives more weight to the
# geometric mean
approx1 <- 1 / 3 * (x + y) / 2 + 2 / 3 * sqrt(x * y)
approx2 <- ((x + y) / 2)^(1 / 3) * (sqrt(x * y))^(2 / 3)
approx1 - logmean(x, y) # always a positive approximation error
#> [1] 1.932965e-02 3.260257e-03 3.420021e-04 3.852275e-06
approx2 - logmean(x, y) # a negative approximation error
#> [1] -6.436118e-02 -1.212079e-02 -1.336412e-03 -1.537117e-05
# A better approximation
correction <- (log(x / y) / pi)^4 / 32
approx1 / (1 + correction) - logmean(x, y)
#> [1] -8.576372e-04 1.064231e-04 2.148563e-05 2.877344e-07
#---- Some identities ----
# A useful identity for turning an additive change into a proportionate
# change
all.equal(logmean(x, y) * log(x / y), x - y)
#> [1] TRUE
# Works for other orders, too
r <- 2
all.equal(
generalized_logmean(r)(x, y)^(r - 1) * (r * (x - y)),
(x^r - y^r)
)
#> [1] TRUE
# Some other identities
all.equal(
generalized_logmean(-2)(1, 2),
(harmonic_mean(1:2) * geometric_mean(1:2)^2)^(1 / 3)
)
#> [1] TRUE
all.equal(
generalized_logmean(0.5)(1, 2),
(arithmetic_mean(1:2) + geometric_mean(1:2)) / 2
)
#> [1] TRUE
all.equal(
logmean(1, 2),
geometric_mean(1:2)^2 * logmean(1, 1 / 2)
)
#> [1] TRUE
#---- Integral representations of the logarithmic mean ----
logmean(2, 3)
#> [1] 2.466303
integrate(function(t) 2^(1 - t) * 3^t, 0, 1)$value
#> [1] 2.466303
1 / integrate(function(t) 1 / (2 * (1 - t) + 3 * t), 0, 1)$value
#> [1] 2.466303