Generalized Linear Model
The maximum likelihood estimate can also be obtained as a generalized
linear model. vmgeom_glm() fits the same likelihood through
glm(), which becomes useful as soon as the r-value is no
longer assumed to be constant.
The response is passed as a two-column matrix holding the meteor
magnitude and the corresponding limiting magnitude. Supplying the
limiting magnitude as part of the response ensures that it is subsetted
together with the remaining data whenever glm() drops rows.
The observed frequencies enter as weights.
# rows with a frequency of zero carry no information
magnitudes_obs <- subset(magnitudes, magnitudes$Freq > 0)
result_glm <- vmgeom_glm(
cbind(magn, lim_magn) ~ 1,
data = magnitudes_obs,
weights = magnitudes_obs$Freq
)
The link function is \(\eta =
\mathrm{logit}(1/r) = -\log(r-1)\). The quantity \(1/r\) is the factor by which the observable
rate drops when the limiting magnitude is lowered by one, and therefore
lies between \(0\) and \(1\), which makes a logit link the natural
choice. It also enforces \(r > 1\)
for every finite \(\eta\), a constraint
the geometric model requires.
Rather than back-transforming the coefficients by hand,
predict() returns the population index directly. Both
methods agree, which is expected: without covariates this is an
intercept-only fit of the same likelihood.
# the model has no covariates, so any single row predicts the global r-value
newdata <- data.frame(row.names = "")
print(c(
"ML" = result_ml$par,
"GLM" = as.numeric(predict(result_glm, newdata, type = "r"))
))
#> ML GLM
#> 2.344528 2.344528
The same holds for the uncertainty. The maximum likelihood standard
deviation comes from the curvature of the log-likelihood, the one of the
linear model from the standard error of \(\eta\), transformed with the delta
method:
print(c(
"ML" = sqrt(1 / result_ml$hessian[1][1]),
"GLM" = as.numeric(predict(result_glm, newdata, type = "r", se.fit = TRUE)$se.fit)
))
#> ML GLM
#> 0.1310968 0.1310969
For a single r-value the likelihood approach is the cheaper one. On
the data set used here it is roughly twice as fast as the linear model,
since glm() carries the overhead of its iterative fitting
on top of the same likelihood. Both are a matter of milliseconds, and
for aggregated frequencies the cost is governed by the number of
magnitude classes rather than by the number of meteors.
Predictions are available on several scales — "r",
"inv_r", "log_r" and "link" — and
se.fit = TRUE adds confidence limits in addition to the
standard error. The limits are formed on the link scale and then
transformed, so they never fall below 1.
print(predict(result_glm, newdata, se.fit = TRUE))
#> $fit
#>
#> 2.344528
#>
#> $se.fit
#>
#> 0.1310969
#>
#> $lwr
#>
#> 2.110642
#>
#> $upr
#>
#> 2.627666
The real benefit appears once covariates are added. Because the fit
is an ordinary glm object, summary(),
anova(), AIC()/BIC() and
select_knots() can all be used. A model in which the
r-value varies with solar longitude would be written as
vmgeom_glm(cbind(magn, lim_magn) ~ sl, data = magnitudes, weights = Freq)
and a smooth activity-dependent profile as
vmgeom_glm(cbind(magn, lim_magn) ~ splines::ns(sl, df = 3), data = magnitudes, weights = Freq)
Since the family reports the exact log-likelihood of the fitted
distribution, AIC() also allows a comparison against a
different magnitude distribution altogether. Fitting the ideal
distribution to the same data with vmideal_glm() gives:
result_ideal <- vmideal_glm(
cbind(magn, lim_magn) ~ 1,
data = magnitudes_obs,
weights = magnitudes_obs$Freq
)
print(c("geometric" = AIC(result_glm), "ideal" = AIC(result_ideal)))
#> geometric ideal
#> 353.1941 350.6568
Each model has one parameter, so the comparison reduces to their
likelihoods. The ideal distribution attains the lower value here, though
by a margin that this excerpt of 97 meteors does not resolve — a
difference of a few units is weak evidence, and the two models are hard
to tell apart at the limiting magnitudes involved. Such a comparison is
only valid when both fits use the same observations and the same
response.
Covariates can also be accommodated on a different basis, by
transforming the magnitudes first; see the section on the
variance-stabilizing transformation below.
By Rate
This method asks how strongly the observed rate would drop if the
limiting magnitude were lowered by one. Under the model of the
introduction, the expected number of meteors recorded at limiting
magnitude \(m_{\mathrm{lim}}\) is
proportional to
\[
S(m_{\mathrm{lim}}) \;=\; \sum_{m \in \mathbb{Z}} f(m_{\mathrm{lim}} -
m)\, r^{m},
\]
each magnitude class contributing its population share \(r^{m}\) weighted by the probability \(f(\cdot)\) of actually perceiving it.
Lowering the limiting magnitude by one and forming the ratio gives the
relative drop of the rate,
\[
\frac{S(m_{\mathrm{lim}} - 1)}{S(m_{\mathrm{lim}})} \;=\; \frac{1}{r},
\]
which turns out to be exactly \(1/r\) — independent of \(m_{\mathrm{lim}}\) and of the shape of
\(f(\cdot)\). Measuring this drop is
what gives the method its name.
The rates themselves are not observed, but the same ratio can be
formed per meteor. For a meteor of magnitude \(m\) observed at limiting magnitude \(m_{\mathrm{lim}}\) let
\[
a(m, m_{\mathrm{lim}}) \;=\; \frac{f(m_{\mathrm{lim}} - m -
1)}{f(m_{\mathrm{lim}} - m)},
\]
the factor by which its chance of being perceived would shrink if the
limiting magnitude were lowered by one, and write \(a_i = a(m_i, m_{\mathrm{lim},i})\) for the
\(i\)-th observed meteor. This leads
to
\[
\mathbb{E}\!\left[\frac{1}{n}\sum_{i=1}^{n} a_i\right] = \frac{1}{r},
\]
so averaging \(a\) over the observed
meteors estimates \(1/r\) without bias.
The derivation below proves both identities.
Derivation
Setting. Fix a limiting magnitude \(m_{\mathrm{lim}}\) and let \(M\) denote the magnitude of a randomly
selected observed meteor. Its distribution is
\[
P[M = m] \;=\; \frac{f(m_{\mathrm{lim}} - m)\,
r^{m}}{S(m_{\mathrm{lim}})}.
\]
Two properties of the perception probability are used, and nothing
else:
- \(f(x) = 0\) for \(x \le -0.5\), that is, a meteor fainter
than the limiting magnitude is never perceived. Hence \(f(m_{\mathrm{lim}} - m) = 0\) for all \(m \ge m_{\mathrm{lim}} +
0.5\), and the sum is bounded above.
- \(0 \le f \le 1\). Together with
\(r > 1\) this bounds the sum from
below as well, since \(\sum_{m \le 0}
r^{m}\) converges. So \(0 <
S(m_{\mathrm{lim}}) < \infty\) and the distribution is
proper.
Let \(m_{\max}\) be the largest
class with \(f(m_{\mathrm{lim}} - m_{\max})
> 0\). The support of \(M\)
is contained in \(\{m \le m_{\max}\}\),
and on it the denominator of \(a(m,
m_{\mathrm{lim}})\) is strictly positive, so \(a(M, m_{\mathrm{lim}})\) is well defined
almost surely.
Claim. \(\mathbb{E}\bigl[a(M, m_{\mathrm{lim}})\bigr] =
1/r\) for every \(m_{\mathrm{lim}}\).
Proof. Writing out the expectation and cancelling
the perception probability, which is legitimate because it is non-zero
on the support,
\[
\begin{aligned}
\mathbb{E}\bigl[a(M, m_{\mathrm{lim}})\bigr]
&= \sum_{m \le m_{\max}} \frac{f(m_{\mathrm{lim}} - m)\,
r^{m}}{S(m_{\mathrm{lim}})}
\cdot \frac{f(m_{\mathrm{lim}} - m - 1)}{f(m_{\mathrm{lim}} - m)}
\\[5pt]
&= \frac{1}{S(m_{\mathrm{lim}})} \sum_{m \le m_{\max}}
f(m_{\mathrm{lim}} - m - 1)\, r^{m}.
\end{aligned}
\]
The remaining sum is by definition \(S(m_{\mathrm{lim}} - 1)\), so that \(\mathbb{E}\bigl[a(M, m_{\mathrm{lim}})\bigr] =
S(m_{\mathrm{lim}} - 1) / S(m_{\mathrm{lim}})\): the expectation
of \(a\) is precisely the relative drop
of the rate stated above.
Before shifting the index, note that its uppermost term already
vanishes. At \(m = m_{\max}\) the
numerator is evaluated one class beyond the faintest perceptible
one,
\[
f\bigl(m_{\mathrm{lim}} - m_{\max} - 1\bigr) = 0,
\]
so the summation effectively stops at \(m_{\max} - 1\). This is the step where the
vanishing perception probability enters, and it is what makes the shift
below land exactly on the range of \(S(m_{\mathrm{lim}})\). Replacing \(m\) by \(m -
1\) therefore gives
\[
\begin{aligned}
\sum_{m \,\le\, m_{\max}} f(m_{\mathrm{lim}} - m - 1)\, r^{m}
&= \sum_{m \,\le\, m_{\max} - 1} f(m_{\mathrm{lim}} - m - 1)\, r^{m}
\\[5pt]
&= \sum_{m \,\le\, m_{\max}} f(m_{\mathrm{lim}} - m)\, r^{m - 1}
\\[5pt]
&= \frac{1}{r} \sum_{m \,\le\, m_{\max}} f(m_{\mathrm{lim}} - m)\,
r^{m} \\[5pt]
&= \frac{S(m_{\mathrm{lim}})}{r},
\end{aligned}
\]
where the last sum is \(S(m_{\mathrm{lim}})\) because \(f(m_{\mathrm{lim}} - m) = 0\) for every
\(m > m_{\max}\), so extending it to
all integers adds nothing. Both series converge absolutely, since \(r^{m} \to 0\) as \(m \to -\infty\), which is what permits the
shift. Altogether
\[
\begin{aligned}
\mathbb{E}\bigl[a(M, m_{\mathrm{lim}})\bigr]
&= \frac{1}{S(m_{\mathrm{lim}})} \cdot \frac{S(m_{\mathrm{lim}})}{r}
\\[5pt]
&= \frac{1}{r}. \qquad \blacksquare
\end{aligned}
\]
Along the way this also proves the rate identity claimed above, since
the intermediate result was \(S(m_{\mathrm{lim}} - 1) =
S(m_{\mathrm{lim}})/r\).
Finally, the second identity follows. For \(n\) meteors observed under limiting
magnitudes \(m_{\mathrm{lim},1}, \dots,
m_{\mathrm{lim},n}\), each \(a_i\) satisfies \(\mathbb{E}[a_i] = 1/r\) by the claim, since
it holds for every limiting magnitude separately. Linearity of
the expectation then gives
\[
\mathbb{E}\!\left[\frac{1}{n}\sum_{i=1}^{n} a_i\right] = \frac{1}{r}.
\]
Properties of the estimator
The sample mean of \(a\) is thus an
unbiased estimator of \(1/r\) —
exactly, not merely asymptotically, and irrespective of whether the
observations share the same limiting magnitude. This is what allows
observations made under different conditions to be pooled. The estimator
is inexpensive to compute and directly comparable to the
likelihood-based approaches discussed earlier.
Unbiasedness for \(1/r\) does not
carry over to \(r\) itself, however:
\(1/\bar{a}\) is a non-linear function
of \(\bar{a}\), which is why the delta
method is applied below.
A second property concerns the bright end. Once a meteor is bright
enough that both perception probabilities are close to 1.0,
the ratio approaches 1 as well, and it no longer matters
how bright the meteor actually was. At a limiting magnitude of
6.0 the ratio is already 0.999 for a magnitude
\(-3\) meteor and indistinguishable
from 1 for anything brighter. Such meteors are therefore
effectively pooled into a single class, and the information about their
exact brightness is lost.
In practice this is not a drawback but an advantage. Meteors that
bright are rare — at \(r = 2.5\) and
the same limiting magnitude, magnitude \(-3\) and brighter accounts for less than
one percent of all observed meteors — so little information is
discarded. In exchange, the pooling keeps a single very bright meteor
from acting as an outlier, precisely in the part of the magnitude range
where one observation would otherwise carry a lot of weight. The same
applies to the variance-stabilizing transformation below, which is built
on this ratio.
Estimation
We estimate \(r\) as \(1/\bar{a}\) from the sample mean \(\bar{a}\). Since \(r\) is a non-linear function of \(\bar{a}\), we apply the delta method to
correct for bias and to compute the variance of \(r\). The delta method uses a Taylor
expansion around the mean of \(a\) to
approximate the distribution of \(r\).
Note that the variance of \(a\)
entering it depends on \(r\) itself;
the section on the variance-stabilizing transformation below starts from
this same statistic and removes that dependency.
result_rate <- with(magnitudes, {
N <- sum(Freq)
a <- vmperception(lim_magn - magn - 1) / vmperception(lim_magn - magn)
a_mean <- as.numeric(weighted.mean(a, w = Freq))
a_var <- as.numeric(cov.wt(cbind(a), wt = Freq)$cov) / N
# apply the delta method and return the result
list(
"mean" = 1 / a_mean - a_var / a_mean^3,
"var" = a_var / a_mean^4
)
})
This gives the expected value and the variance of the r-value:
print(result_rate$mean) # mean of r
#> [1] 2.323974
print(result_rate$var) # variance of r
#> [1] 0.01273307
Using the bootstrap method, it can be assessed whether the mean is
normally distributed.
# Bootstrapping Method
r_means <- with(magnitudes, {
N <- sum(Freq)
a <- vmperception(lim_magn - magn - 1) / vmperception(lim_magn - magn)
replicate(50000, {
s <- sample(a, size = N, replace = TRUE, prob = Freq)
1 / mean(s)
})
})
The graphical representation indicates that this is indeed
approximately the case.
with(new.env(), {
r_sd <- sqrt(result_rate$var)
r_min <- result_rate$mean - 3 * r_sd
r_max <- result_rate$mean + 3 * r_sd
r <- subset(r_means, r_means > r_min & r_means < r_max)
brks <- seq(min(r) - 0.02, max(r) + 0.02, by = 0.02)
hist(r,
breaks = brks,
col = "skyblue",
border = "black",
main = "Histogram of mean r",
xlab = "r",
xaxt = "n",
ylab = "count"
)
xlabels <- seq(min(round(r, 1)) - 0.1, max(round(r, 1)) + 0.1, by = 0.1)
axis(
side = 1,
at = xlabels,
labels = sprintf("%.1f", xlabels)
)
abline(v = result_rate$mean, col = "red", lwd = 1)
})
Variance-Stabilizing Transformation as an Alternative Method
This method is a refinement of the rate-based one. It uses the same
quantity \(a(m, m_{\mathrm{lim}})\),
but rescales it so that its variance no longer depends on the parameter
r. It is likewise far cheaper than estimation based on the
maximum likelihood principle.
Recall what was established in the previous section. The sample mean
of \(a\) estimates \(1/r\) without bias, exactly and for pooled
limiting magnitudes alike, which is a strong property. Its weakness is
the second moment: the variance of \(a\) is itself a function of \(r\). At a limiting magnitude of
6.0 it falls from about 0.066 at \(r = 1.4\) to 0.033 at \(r = 4\), a factor of two across the range
of interest. Two consequences follow. The precision of the estimate
cannot be judged without already knowing \(r\), and the delta method is needed to
convert \(\bar{a}\) into an estimate of
\(r\) at all.
The variance-stabilizing transformation removes exactly this
weakness. It maps \(a\) onto a scale on
which the variance is approximately 1.0 regardless of \(r\) — over \(1.4
\le r \le 4\) it stays within about 0.16 of that
value, comfortably covering the \(1.7 \le r
\le 3.3\) met in visual meteor work. Outside that window the
stabilization degrades on both sides, which bounds the use of the
transformed magnitudes as a response in linear models, but not the
reading of a single \(r\) off their
mean. vmgeom_vst_from_magn() performs the mapping and
vmgeom_vst_to_r() maps the mean of the transformed
magnitudes back onto \(r\).
The mapping is monotonically decreasing: a bright meteor is far from
the limiting magnitude and receives a large tm. Towards the
bright end it runs into a ceiling, since \(a\) cannot exceed 1 and the
transformed value is therefore bounded by what the transformation makes
of that limit. At the faint end it reaches zero half a magnitude below
the limiting magnitude, where a meteor is no longer perceived. Shown
here for a limiting magnitude of 6.0.
with(new.env(), {
lim_magn <- 6.0
m <- seq(-2, 5.5, 0.05)
plot(m, vmgeom_vst_from_magn(m, lim_magn),
type = "l",
col = "blue",
xlab = "m",
ylab = "tm"
)
})
The way back needs no limiting magnitude at all, since the statistic
it rests on does not depend on one. On logarithmic axes the relation is
nearly a straight line, which is the additivity discussed further
below.
with(new.env(), {
tm <- seq(1.5, 4.3, 0.01)
plot(tm, vmgeom_vst_to_r(tm),
type = "l",
col = "blue",
log = "xy",
xlab = "tm",
ylab = "r",
xaxt = "n",
yaxt = "n"
)
axis(1, at = c(1.5, 2, 3, 4), labels = c("1.5", "2", "3", "4"))
axis(2, at = c(1, 1.5, 2, 3, 4), labels = c("1", "1.5", "2", "3", "4"))
})
The price is paid on the first moment. The transformation is
non-linear in \(a\), and the mean of a
non-linear function is not that function of the mean, so the exact
unbiasedness of the rate-based estimator is lost. This deviation is
systematic: it does not average out, and it does not shrink as the
sample grows.
Its cause is a slight dependence of the transformed mean on the
limiting magnitude, which the statistic \(a\) itself does not have. Over \(1.7 \le r \le 3.3\), the range met in
visual meteor work, it amounts to at most 0.8% of \(r\).
What matters in practice is how that compares with the random error,
and this depends on the number of meteors. The random error shrinks as
the sample grows while the systematic one does not, so the two meet at a
certain sample size: for a few hundred meteors the random error still
dominates by a factor of five or more, at a few thousand by roughly two,
and only beyond some ten thousand does the systematic part set the
accuracy — no matter how many more meteors are added. Visual meteor data
sets of that size are the exception, so for most analyses the deviation
can be ignored.
Where it does matter — for a very large data set, or wherever an
unbiased estimate counts for more than a constant variance — the
rate-based method of the previous section remains exactly unbiased for
\(1/r\), and vmgeom_glm()
fits the exact likelihood.
Everything the rate-based method offers is retained. Since the
transformation is a monotone rescaling of \(a\), the pooling of bright meteors
described above carries over unchanged. And because \(a\) lies between 0 and
1, the transformed values are bounded as well, with the
upper bound corresponding to \(r = 1\);
the estimate can consequently never fall below the value the geometric
model requires.
The same monotonicity matters when the transformation is used inside
a predictive model. vmgeom_vst_to_r() is defined on the
whole range of values the transformation can produce and returns
NA only for values it cannot. Should sparse data drive a
prediction into an implausible region, the result is therefore a value
that is too large but still ordered, rather than a missing one — usually
far easier to handle downstream than an NA.
The delta method is applied below in this section as well, since the
back-transformation is non-linear here too. What changes is the variance
that enters it: a constant one instead of a value that depends on the
very parameter being estimated.
Beyond stabilizing the variance, the transformation also puts the
parameter on a scale with a simple structure. Classically,
transformations of this kind are judged by three criteria, which Box and
Cox (1964, An Analysis of Transformations, J. R. Stat. Soc. B
26, 211–252) treat together: constant variance, an approximately
symmetric distribution, and additivity, meaning that effects
combine linearly on the transformed scale. The third one is what the
transformation is built around. It belongs to the Box-Cox family itself,
applied to \(a\), and \(a\) estimates \(1/r\), so
\[
\log t \;=\; \mathrm{const.} - \beta \, \log r
\]
with \(\beta > 0\): the relation
is close to linear on the logarithmic scale. A multiplicative change of
the population index therefore becomes an approximately additive shift
of \(\log t\). Over the range of \(r\) the calibration covers this is what the
mean of the transformed magnitudes follows, to the accuracy stated
above, and it is what the calibration of vmgeom_vst_to_r()
rests on. The departure from exact linearity is slight but systematic,
and the back-transformation accounts for it explicitly — see below.
This is what makes the transformed magnitudes usable as a response in
ordinary linear models, with residuals that are homoscedastic by
construction — precisely what such models assume.
vmgeom_vst_lm() fits them, and the estimation below is
already of that form: taking the mean of the transformed magnitudes
is an intercept-only linear model. A profile in which \(r\) varies with a covariate only extends
the right-hand side of the formula, so that
vmgeom_vst_lm(cbind(magn, lim_magn) ~ sl, data = magnitudes, weights = Freq)
fits a solar-longitude dependence in the same way
vmgeom_glm() does above. The estimated coefficients
translate back into changes of \(\log
r\) up to the factor \(\beta\).
The back-transformation is \(\log r = c + d
\log t + e\,t\), so that factor is \(\beta = -1/(d + e\,t)\), which over the
range met in practice runs from about 0.73 to
0.81. The term in \(t\) is
what captures the curvature of the relation; a straight line through it
would leave a systematic error about 1.7 times as large,
and it is the reason vmgeom_vst_to_r() reproduces \(r = 1\) exactly at the upper bound of
tm.
That makes this approach a lightweight alternative to
vmgeom_glm(), and the choice between them is a trade-off.
In favour of the transformation is its cost: on the data set used here a
single estimate takes well under a millisecond, some two orders of
magnitude less than the generalized linear model, which iterates over
the full likelihood. Since the transformation is applied once per meteor
and everything afterwards is ordinary least squares, this advantage
grows with the size of the data set, and the whole toolbox built around
lm() becomes available.
Against it stands a loss of accuracy, though not where one would
expect it. The dispersion of the estimate is practically that of the
maximum likelihood estimate — in a simulation with n = 500
the two standard deviations agree to within a few percent over \(2 \le r \le 4\). What the transformation
costs is not variance but the systematic deviation of the
back-transformation described above, which adds to the total error and,
unlike the random part, does not shrink with the sample size. At
n = 500 and \(r = 3\) it
raises the root mean squared error by roughly 3% over that
of the maximum likelihood estimate; for larger samples its share grows.
vmgeom_glm(), fitting the exact likelihood, carries no such
deviation and additionally provides the inferential machinery of a
glm object — summary(), anova(),
information criteria. The transformation is therefore preferable for
quick estimates and exploratory work, the generalized linear model
whenever accuracy matters — including, unlike the usual trade-off, for
large data sets, where the systematic part comes to dominate.
A constant variance carries a few further consequences worth
noting:
- Standard errors and confidence intervals no longer depend on the
unknown r-value, so they are equally well calibrated across the entire
parameter space,
- Since the sample mean is homoscedastic, the classical results of
linear estimation theory (BLUE properties, Gauss–Markov theorem) apply
without further transformation,
- Test statistics such as likelihood ratio or score tests follow a
more uniform distribution, because the variance need not be treated as
an additional unknown. This improves test power and simplifies
asymptotic approximations.
The resulting procedure is straightforward: all that is needed is the
mean of the transformed meteor magnitudes, from which an estimate of the
parameter r is obtained. vmgeom_vst_lm() does
both in one step. It applies the transformation and fits the result as
the intercept-only linear model announced above, so that the step to a
covariate is a matter of extending the formula. The response is passed
as cbind(magn, lim_magn) and the observed frequencies as
weights, exactly as for vmgeom_glm().
result_vs <- vmgeom_vst_lm(
cbind(magn, lim_magn) ~ 1,
data = magnitudes,
weights = magnitudes$Freq
)
predict() converts the fitted mean back to the
r scale. It applies the delta method, which accounts for
the nonlinearity of the back-transformation, and returns the standard
error alongside the estimate.
# the model has no covariates, so any single row predicts the global r-value
newdata <- data.frame(row.names = "")
prediction <- predict(result_vs, newdata, se.fit = TRUE, bias_correction = TRUE)
Two of the arguments deserve a comment. Passing the frequencies as
weights is the correct and complete way to supply
aggregated counts — nothing is lost by it, and there is no separate
mechanism one ought to be using instead. The coefficients it produces
are exactly those of a data set holding one row per meteor.
What differs is only how the residual variance is read back.
R has no notion of frequency weights in lm(): weights are
taken to express precision, so the residual variance is normalized by
the number of rows. summary() therefore reports a standard
error that refers to the 48 magnitude classes rather than to the 97
meteors behind them — here about 40% too large, with the
confidence interval correspondingly too wide. The same applies to
glm(family = gaussian()).
predict() takes care of this: it derives the residual
scale from the sum of the weights and thus reports the standard error of
the meteors. Nothing needs to be done differently when calling the
function — but the standard errors of summary(result_vs)
should not be used, since they carry the row-based normalization.
Expanding the rows into one meteor each would settle the question at
the source, and for integer counts it gives exactly the same result as
the correction applied here. Visual observations, however, may contain
half meteors, which makes that route unavailable.
bias_correction = TRUE adds the second-order term of the
delta method. It is not the default — predict() methods in
R report a first-order standard error, and so does
vmgeom_glm() above — but it is what this section is about.
The term is small here, well below the random error at this sample size,
and it corrects the curvature of the back-transformation, not
the systematic deviation described further above. That deviation is a
property of the calibration and does not shrink as the sample grows.
Thus, one obtains the mean and the variance of the mean of
r.
result_vs <- list(mean = unname(prediction$fit), var = unname(prediction$se.fit)^2)
print(paste("r mean:", result_vs$mean))
#> [1] "r mean: 2.26269452053892"
print(paste("r var:", result_vs$var))
#> [1] "r var: 0.0103330446553896"
Using the bootstrap method, it can be assessed whether the mean is
normally distributed.
# Bootstrapping Method
r_means <- with(magnitudes, {
N <- sum(Freq)
tm <- vmgeom_vst_from_magn(magn, lim_magn)
replicate(50000, {
s <- sample(tm, size = N, replace = TRUE, prob = Freq)
vmgeom_vst_to_r(mean(s))
})
})
The graphical representation indicates that this is indeed
approximately the case.
with(new.env(), {
r_sd <- sqrt(result_vs$var)
r_min <- as.vector(result_vs$mean - 3 * r_sd)
r_max <- as.vector(result_vs$mean + 3 * r_sd)
r <- subset(r_means, r_means > r_min & r_means < r_max)
brks <- seq(min(r) - 0.02, max(r) + 0.02, by = 0.02)
hist(r,
breaks = brks,
col = "skyblue",
border = "black",
main = "Histogram of mean r",
xlab = "r",
xaxt = "n",
ylab = "count"
)
xlabels <- seq(min(round(r, 1)) - 0.1, max(round(r, 1)) + 0.1, by = 0.1)
axis(
side = 1,
at = xlabels,
labels = sprintf("%.1f", xlabels)
)
abline(v = result_vs$mean, col = "red", lwd = 1)
})
Residual Analysis
So far, we have operated under the assumption that the real
distribution of meteor magnitudes is exponential and that the perception
probabilities are accurate. We now use the Chi-Square goodness-of-fit
test to check whether the observed frequencies match the expected
frequencies. Then, using the estimated r-value, we retrieve the relative
frequencies p for each observation and add them to the data
frame magnitudes:
magnitudes$p <- with(magnitudes, dvmgeom(m = magn, lm = lim_magn, result_rate$mean))
We must also consider the probabilities for the magnitude class with
the brightest meteors.
magn_min <- min(magnitudes$magn)
The smallest magnitude class magn_min is -6. In
calculating the probabilities, we assume that the magnitude class -6
contains meteors that are either brighter or equally bright as -6 and
thus use the function pvmgeom() to determine their
probability.
idx <- magnitudes$magn == magn_min
magnitudes$p[idx] <- with(
magnitudes[idx, ],
pvmgeom(m = magn + 1L, lm = lim_magn, result_rate$mean, lower.tail = TRUE)
)
This ensures that the probability of observing a meteor of any given
magnitude is 100%. This is known as the normalization condition.
Accordingly, the Chi-Square goodness-of-fit test will fail if this
condition is not met.
We now create the contingency table magnitutes_observed
for the observed meteor magnitudes and its margin table.
magnitutes_observed <- xtabs(Freq ~ magn_id + magn, data = magnitudes)
magnitutes_observed_mt <- margin.table(magnitutes_observed, margin = 2)
print(magnitutes_observed_mt)
#> magn
#> -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
#> 0.0 0.0 0.0 0.0 3.0 4.0 7.0 10.0 23.0 26.5 20.0 3.0 0.5
Next, we check which magnitude classes need to be aggregated so that
each contains at least 10 meteors, allowing us to perform a Chi-Square
goodness-of-fit test.
The last output shows that meteors of magnitude class 0
or brighter must be combined into a magnitude class 0-.
Meteors with a brightness less than 4 are grouped here in
the magnitude class 4+, and a new contingency table
magnitudes.observed is created:
magnitudes$magn[magnitudes$magn <= 0] <- "0-"
magnitudes$magn[magnitudes$magn >= 4] <- "4+"
magnitutes_observed <- xtabs(Freq ~ magn_id + magn, data = magnitudes)
print(margin.table(magnitutes_observed, margin = 2))
#> magn
#> 0- 1 2 3 4+
#> 14.0 10.0 23.0 26.5 23.5
We now need the corresponding expected relative frequencies
magnitutes_expected <- xtabs(p ~ magn_id + magn, data = magnitudes)
magnitutes_row_freq <- margin.table(magnitutes_observed, margin = 1)
magnitutes_expected <- sweep(magnitutes_expected, 1, magnitutes_row_freq, `*`)
magnitutes_expected <- magnitutes_expected / sum(magnitutes_expected)
print(sum(magnitudes$Freq) * margin.table(magnitutes_expected, margin = 2))
#> magn
#> 0- 1 2 3 4+
#> 15.68943 14.26433 19.69970 21.01861 26.32794
and then carry out the Chi-Square goodness-of-fit test:
chisq_test_result <- chisq.test(
x = margin.table(magnitutes_observed, margin = 2),
p = margin.table(magnitutes_expected, margin = 2)
)
As a result, we obtain the p-value:
chi2_df <- chisq_test_result$parameter - 1
chi2_pval <- pchisq(chisq_test_result$statistic, df = chi2_df, lower.tail = FALSE)
print(chi2_pval)
#> X-squared
#> 0.2906011
If we set the level of significance at 5 percent, then it is clear
that the p-value with 0.2906011 is greater than 0.05. Thus, under the
assumption that the magnitude distribution follows an geometric meteor
magnitude distribution and assuming that the perception probabilities
are correct (i.e., error-free or precisely known), the assumptions
cannot be rejected. However, the converse is not true; the assumptions
may not necessarily be correct. The total count of meteors here is too
small for such a conclusion.
To verify the p-value, we also graphically represent the Pearson
residuals:
chisq_test_residuals <- with(new.env(), {
chisq_test_residuals <- residuals(chisq_test_result)
v <- as.vector(chisq_test_residuals)
names(v) <- names(chisq_test_residuals)
v
})
plot(
chisq_test_residuals,
main = "Residuals of the chi-square goodness-of-fit test",
xlab = "m",
ylab = "Residuals",
ylim = c(-3, 3),
xaxt = "n"
)
abline(h = 0.0, lwd = 2)
axis(1, at = seq_along(chisq_test_residuals), labels = names(chisq_test_residuals))
