--- title: "Parameter Estimation of the Geometric Model of Visual Meteor Magnitudes" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true fig_width: 6 fig_height: 4 vignette: > %\VignetteIndexEntry{Parameter Estimation of the Geometric Model of Visual Meteor Magnitudes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} library(vismeteor) knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Introduction The geometric model of visual meteor magnitudes is a commonly used statistical approach to describe the magnitude distribution of a meteor shower. The observable magnitude distribution of meteors is then $$ P[M = m] \sim \begin{cases} f(m_{\mathrm{lim}} - m)\, r^m, & \text{if } m_{\mathrm{lim}} - m > -0.5,\\[5pt] 0 & \text{otherwise,} \end{cases} $$ where $m_{\mathrm{lim}}$ denotes the limiting (non-integer) magnitude of the observation, and $m$ the integer meteor magnitude. The function $f(\cdot)$ denotes the perception probability function. The estimation of the population index r, briefly called the r-value, is a common task in the evaluation of meteor magnitudes. Here we demonstrate several methods to estimate this parameter. First, we obtain some magnitude observations from the example data set, which also includes the limiting magnitude. ```{r, echo=TRUE, results='hide'} observations <- with(PER_2015_magn$observations, { idx <- !is.na(lim_magn) & sl_start > 135.81 & sl_end < 135.87 data.frame( magn_id = magn_id[idx], lim_magn = lim_magn[idx] ) }) head(observations, 5) # Example values ``` ```{r, echo=FALSE, results='asis'} knitr::kable(head(observations, 5)) ``` Next, the observed meteor magnitudes are matched with the corresponding observations. This is necessary as we need the limiting magnitudes of the observations to determine the r-value. Using ```{r, echo=TRUE, results='hide'} magnitudes <- with(new.env(), { magnitudes <- merge( observations, as.data.frame(PER_2015_magn$magnitudes), by = "magn_id" ) magnitudes$magn <- as.integer(as.character(magnitudes$magn)) subset(magnitudes, (magnitudes$lim_magn - magnitudes$magn) > -0.5) }) head(magnitudes, 5) # Example values ``` we obtain a data frame with the absolute observed frequencies `Freq` for each observation of a magnitude class. The expression `subset(magnitudes, (magnitudes$lim_magn - magnitudes$magn) > -0.5` ensures that meteors fainter than the limiting magnitude are not used if they exist. ```{r, echo=FALSE, results='asis'} knitr::kable(head(magnitudes[magnitudes$Freq > 0, ], 5)) ``` This data frame contains a total of `r sum(magnitudes$Freq)` meteors. This is a sufficiently large number to estimate the r-value. ## Maximum Likelihood Method The maximum likelihood method can be used to estimate the r-value in an asymptotically unbiased manner. For this, the function `dvmgeom()` is needed, which returns the probability density of the observable meteor magnitudes when the r-value and the limiting magnitudes are known. The following algorithm estimates the r-value by maximizing the likelihood with the `optim()` function. The function `ll()` returns the negative log-likelihood, as `optim()` identifies a minimum. ```{r, echo=TRUE, results='hide'} # maximum likelihood estimation (MLE) of r result_ml <- with(magnitudes, { # log likelihood function ll <- function(r) -sum(Freq * dvmgeom(magn, lim_magn, r, log = TRUE)) r_start <- 2.0 # starting value r_lower <- 1.2 # lowest expected value r_upper <- 4.0 # highest expected value # find minimum optim(r_start, ll, method = "Brent", lower = r_lower, upper = r_upper, hessian = TRUE) }) ``` This gives the expected value and the variance of the r-value: ```{r, echo=TRUE} print(result_ml$par) # mean of r print(1 / result_ml$hessian[1][1]) # variance of r ``` We can additionally visualize the likelihood function here. ```{r, echo=TRUE, results='hide'} with(new.env(), { data_plot <- data.frame(r = seq(2.0, 2.8, 0.01)) data_plot$ll <- mapply(function(r) { with(magnitudes, { # log likelihood function sum(Freq * dvmgeom(magn, lim_magn, r, log = TRUE)) }) }, data_plot$r) data_plot$l <- exp(data_plot$ll - max(data_plot$ll)) data_plot$l <- data_plot$l / sum(data_plot$l) brks <- seq(min(data_plot$r) - 0.02, max(data_plot$r) + 0.02, by = 0.02) plot(data_plot$r, data_plot$l, # breaks = brks, type = "l", col = "blue", xlab = "r", xaxt = "n", ylab = "likelihood" ) xlabels <- seq(min(round(data_plot$r, 1)) - 0.1, max(round(data_plot$r, 1)) + 0.1, by = 0.1) axis( side = 1, at = xlabels, labels = sprintf("%.1f", xlabels) ) abline(v = result_ml$par, col = "red", lwd = 1) }) ``` The likelihood function is approximately a normal distribution in this case. This is important in this context because the variance of the estimated r-value is derived from the curvature (the second derivative at the maximum) of the log-likelihood function. ## 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`. ```{r, echo=TRUE, results='hide'} # 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. ```{r, echo=TRUE} # 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")) )) ``` 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: ```{r, echo=TRUE} print(c( "ML" = sqrt(1 / result_ml$hessian[1][1]), "GLM" = as.numeric(predict(result_glm, newdata, type = "r", se.fit = TRUE)$se.fit) )) ``` 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`. ```{r, echo=TRUE} print(predict(result_glm, newdata, se.fit = TRUE)) ``` 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 ```r vmgeom_glm(cbind(magn, lim_magn) ~ sl, data = magnitudes, weights = Freq) ``` and a smooth activity-dependent profile as ```r 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: ```{r, echo=TRUE} 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))) ``` 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 `r sum(magnitudes$Freq)` 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: 1. $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. 2. $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. ```{r, echo=TRUE, results='hide'} 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: ```{r, echo=TRUE} print(result_rate$mean) # mean of r print(result_rate$var) # variance of r ``` Using the bootstrap method, it can be assessed whether the mean is normally distributed. ```{r, echo=TRUE, results='hide'} # 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. ```{r, echo=TRUE, results='show'} 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`. ```{r, echo=TRUE, results='show'} 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. ```{r, echo=TRUE, results='show'} 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 ```r 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()`. ```{r, echo=TRUE, results='show'} 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. ```{r, echo=TRUE, results='show'} # 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 `r nrow(magnitudes[magnitudes$Freq > 0, ])` magnitude classes rather than to the `r sum(magnitudes$Freq)` 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`. ```{r, echo=TRUE, results='show'} result_vs <- list(mean = unname(prediction$fit), var = unname(prediction$se.fit)^2) print(paste("r mean:", result_vs$mean)) print(paste("r var:", result_vs$var)) ``` Using the bootstrap method, it can be assessed whether the mean is normally distributed. ```{r, echo=TRUE, results='hide'} # 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. ```{r, echo=TRUE, results='show'} 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`: ```{r, echo=TRUE, results='hide'} 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. ```{r, echo=TRUE, results='hide'} magn_min <- min(magnitudes$magn) ``` The smallest magnitude class `magn_min` is `r magn_min`. In calculating the probabilities, we assume that the magnitude class `r magn_min` contains meteors that are either brighter or equally bright as `r magn_min` and thus use the function `pvmgeom()` to determine their probability. ```{r, echo=TRUE, results='asis'} 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. ```{r, echo=TRUE} magnitutes_observed <- xtabs(Freq ~ magn_id + magn, data = magnitudes) magnitutes_observed_mt <- margin.table(magnitutes_observed, margin = 2) print(magnitutes_observed_mt) ``` 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: ```{r, echo=TRUE} 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)) ``` We now need the corresponding expected relative frequencies ```{r, echo=TRUE} 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)) ``` and then carry out the Chi-Square goodness-of-fit test: ```{r, echo=TRUE, results='asis'} 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: ```{r, echo=TRUE} chi2_df <- chisq_test_result$parameter - 1 chi2_pval <- pchisq(chisq_test_result$statistic, df = chi2_df, lower.tail = FALSE) print(chi2_pval) ``` If we set the level of significance at 5 percent, then it is clear that the p-value with `r unname(chi2_pval)` 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: ```{r, fig.show='hold'} 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)) ```