--- title: "Parameter Estimation of the Ideal Distribution 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 Ideal Distribution 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 density of the ideal distribution of meteor magnitudes is $$ {\displaystyle f(m) = \frac{\mathrm{d}p}{\mathrm{d}m} = \frac{3}{2} \, \log(r) \sqrt{\frac{r^{3 \, \psi + 2 \, m}}{(r^\psi + r^m)^5}}} $$ where $m$ denotes the continuous (real-valued) meteor magnitude, $r = 10^{0.4} \approx 2.51189 \dots$ is a constant, and $\psi$ is the only parameter of this magnitude distribution. In visual meteor observations, magnitudes are usually estimated as integer values. Hence, this distribution is discrete and its probability mass function is given by $$ P[M = m] \sim \begin{cases} g(m_{\mathrm{lim}} - m) \displaystyle \int\limits_{m-0.5}^{m+0.5} f(u) \, \mathrm{d}u, & \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)$ is the continuous density of the ideal magnitude distribution, and $g(\cdot)$ denotes the perception probability function. Here we demonstrate two methods to estimate the parameter $\psi$. 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 parameter. 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[magnitudes$Freq > 0, ], 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 parameter. ## Maximum Likelihood Method The maximum likelihood method can be used to estimate the parameter in an asymptotically unbiased manner. For this, the function `dvmideal()` is needed, which returns the probability density of the observable meteor magnitudes when the parameter and the limiting magnitudes are known. The following algorithm estimates the parameter 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 psi result_ml <- with(magnitudes, { # log likelihood function ll <- function(psi) -sum(Freq * dvmideal(magn, lim_magn, psi, log = TRUE)) psi_start <- 6.0 # starting value psi_lower <- 4.0 # lowest expected value psi_upper <- 10.0 # highest expected value # find minimum optim(psi_start, ll, method = "Brent", lower = psi_lower, upper = psi_upper, hessian = TRUE) }) ``` This gives the expected value and the variance of the parameter: ```{r, echo=TRUE} psi_mean <- result_ml$par # mean of psi print(psi_mean) psi_var <- 1 / result_ml$hessian[1][1] # variance of psi print(psi_var) ``` `optim()` requires both bounds to be finite, so `psi_upper` cannot be set to infinity. That bound is an assumption, and it is worth testing rather than trusting: an estimate that has merely run into it would be reported like any other. Since `dvmideal()` accepts an infinite $\psi$, the likelihood can be evaluated at the limit itself and compared with the estimate. ```{r, echo=TRUE} ll_at <- function(psi) with(magnitudes, sum(Freq * dvmideal(magn, lim_magn, psi, log = TRUE))) print(c("psi_hat" = ll_at(psi_mean), "Inf" = ll_at(Inf))) ``` The estimate is the better of the two, so the maximum lies within the interval and $\psi$ is identified. Had the limit won instead, the maximum likelihood estimate would be `Inf`: the magnitudes are then geometric with $r = 10^{0.4}$, and no finite $\psi$ describes them better. Raising `psi_upper` would not help in that case, and beyond about `35` it becomes actively harmful here — the likelihood is flat above $\psi \approx 16.5$, and once that plateau dominates the search interval, `optim()` settles on it instead of finding the maximum. We can additionally visualize the likelihood function here. ```{r, echo=TRUE, results='hide'} with(new.env(), { data_plot <- data.frame(psi = seq(4.0, 11, 0.1)) data_plot$ll <- mapply(function(psi) { with(magnitudes, { # log likelihood function sum(Freq * dvmideal(magn, lim_magn, psi, log = TRUE)) }) }, data_plot$psi) data_plot$l <- exp(data_plot$ll - max(data_plot$ll)) data_plot$l <- data_plot$l / sum(data_plot$l) plot(data_plot$psi, data_plot$l, type = "l", col = "blue", xlab = "psi", ylab = "likelihood" ) abline(v = result_ml$par, col = "red", lwd = 1) }) ``` It is clearly visible that the likelihood function is not normally distributed. This distribution even belongs to the class of heavy-tailed distributions. While its maximum is indeed an asymptotically unbiased estimator, this does not hold for the variance. This is important in this context because the variance of the estimated $\psi$-value is derived from the curvature (the second derivative at the maximum) of the log-likelihood function. Therefore, the estimator for the variance of $\psi$ is far too small. ## Variance-Stabilizing Transformation Estimation based on the maximum likelihood principle is computationally demanding. As an alternative to the method of the previous section, a variance-stabilizing transformation can be applied. It maps meteor magnitudes onto a different scale, yielding a distribution whose variance no longer depends on the parameter $\psi$, so that a plain mean of the transformed magnitudes estimates $\psi$. What this buys is a fixed variance of `1.0`, and several things follow from it: - Standard errors can be computed without knowing $\psi$, since the variance is no longer a function of the parameter, - Confidence intervals are equally well calibrated across the entire parameter space, - Observations made at different limiting magnitudes become directly comparable, because they now live on a common scale with the same dispersion, - Test statistics exhibit a more uniform distribution, since the variance does not need to be treated as an additional unknown. The mean of the transformed magnitudes is therefore a well-conditioned estimator, and its dispersion no longer depends on the parameter being estimated. It is worth being precise about what this does and does not mean: the sample mean becomes easy to work with, but the transformed magnitudes do not thereby become a quantity that is linear in covariates. The next paragraphs return to this point, because the distinction decides which tool applies. The resulting procedure is straightforward: it suffices to compute the mean of the transformed meteor magnitudes, from which an estimate of the parameter $\psi$ is obtained. `vmideal_vst_from_magn()` performs the mapping. It is monotonically decreasing: a bright meteor is far from the limiting magnitude and receives a large `tm`, a faint one a small value. Shown here for a limiting magnitude of `6.0`. ```{r, echo=TRUE, results='show'} with(new.env(), { lim_magn <- 6.0 m <- seq(-4, 6, 0.05) plot(m, vmideal_vst_from_magn(m, lim_magn), type = "l", col = "blue", xlab = "m", ylab = "tm" ) abline(h = 0.0, lwd = 1) }) ``` Note that `tm` becomes negative for meteors close to the limiting magnitude. This is where the transformation runs out: `vmideal_vst_to_psi()` maps the mean back onto $\psi$, and it is defined from `tm = 0.02` upwards, reporting `Inf` below that value. The two directions therefore cover different ranges of `tm`, which is not a defect but the statement that faint meteors no longer resolve a finite $\psi$. ```{r, echo=TRUE, results='show'} with(new.env(), { lim_magn <- 6.0 tm <- seq(0.02, 8.22, 0.02) plot(tm, vmideal_vst_to_psi(tm, lim_magn), type = "l", col = "blue", xlab = "tm", ylab = "psi" ) abline(h = lim_magn, col = "red", lwd = 1) }) ``` The red line marks the limiting magnitude itself. Note that $\psi$ decreases as `tm` grows: the region above the line, where $\psi$ exceeds the limiting magnitude and the mean of the distribution begins to flatten, is therefore the narrow band at small `tm` on the left. That is where the estimate loses its resolution, and it is also why a larger `tm` corresponds to a smaller $\psi$ — a relation the confidence bound below relies on. Note that the variance-stabilizing transformation yields a mean that can be directly analyzed, but converting it back to the $\psi$ parameter requires applying the delta method. This accounts for the nonlinearity of the transformation and provides appropriate uncertainty estimates for $\psi$. Because the back-transformation returns $\psi$ relative to a limiting magnitude, a reference point has to be fixed before anything can be estimated. The transformed magnitude `tm` measures how far $\psi$ lies behind the limiting magnitude of the observation, so a value of `tm` on its own does not identify $\psi$ — only together with the limiting magnitude it refers to. The mean limiting magnitude of the data set is the natural choice, being the point at which the estimate is best supported. ```{r, echo=TRUE} lim_magn_mean <- with(magnitudes, { N <- sum(Freq) sum(Freq * lim_magn) / N }) print(paste("lim_magn_mean:", lim_magn_mean)) ``` ```{r, echo=TRUE, results='show'} tm_mean <- with(magnitudes, { N <- sum(Freq) tm <- vmideal_vst_from_magn(magn, lim_magn) tm_mean <- sum(Freq * tm) / N tm_var <- sum(Freq * (tm - tm_mean)^2) / (N - 1) tm_mean_var <- tm_var / N list("val" = tm_mean, "var" = tm_mean_var, "sd" = sqrt(tm_mean_var)) }) ``` Thus, one obtains the mean and the variance of the mean of `tm`. ```{r, echo=TRUE, results='show'} print(paste("tm mean:", tm_mean$val)) print(paste("tm var:", tm_mean$var)) ``` Using the bootstrap method, it can be assessed whether the mean is normally distributed. ```{r, echo=TRUE, results='show'} tm_means <- with(magnitudes, { N <- sum(Freq) tm <- vmideal_vst_from_magn(magn, lim_magn) replicate(50000, { mean(sample(tm, size = N, replace = TRUE, prob = Freq)) }) }) ``` The graphical representation indicates that this is indeed approximately the case. ```{r, echo=TRUE, results='show'} with(new.env(), { tm_min <- tm_mean$val - 3 * tm_mean$sd tm_max <- tm_mean$val + 3 * tm_mean$sd tm_means <- subset(tm_means, tm_means > tm_min & tm_means < tm_max) brks <- seq(min(tm_means) - 0.02, max(tm_means) + 0.02, by = 0.02) hist(tm_means, breaks = brks, col = "skyblue", border = "black", main = "Histogram of mean tm", xlab = "tm", ylab = "count", xaxt = "n" ) axis(1, at = seq(round(min(brks), 1), round(max(brks), 1) + 0.1, by = 0.1)) abline(v = 0, col = "red", lwd = 1) }) ``` A mean value approaching `0.0` implies that $\psi$ lies at infinity, and `vmideal_vst_to_psi()` reports `Inf` once it has fallen below `0.02`, where no finite $\psi$ is resolved any more. Negative values can be interpreted as a kind of “beyond infinity”, which is not meaningful. There are two possible explanations: 1. The distribution is not ideal, i.e., the observations do not fit the model as described above. 2. Random variation led to this result. This data set is a case in point. The mean lies only about one standard deviation above `0.0`, and roughly a quarter of the bootstrap replicates fall below the threshold at which no finite $\psi$ is resolved. The point estimate is therefore better read as a median than as an expected value, and the confidence bound given below is the more honest statement. ```{r, echo=TRUE} print(paste("mean psi:", vmideal_vst_to_psi(tm_mean$val, lim_magn_mean))) ``` No second-order term of the delta method is applied here, and this is the same circumstance seen from another side. `vmideal_vst_to_psi()` supplies the necessary derivatives through `deriv_degree`, as `vignette("vmgeom")` uses them, but such a correction presupposes that the back-transformation is nearly straight over the range the estimate covers. Here it is not: the second-order term amounts to some nine percent of $\psi$, and the lower end of the interval falls below the threshold where $\psi$ ceases to be finite. A correction of that size no longer corrects anything — it reports that the expansion it comes from does not apply. Where the estimate is well separated from the threshold the term is legitimate, and then it is too small to matter. This is the case the transformation is made for: a single, global $\psi$ that does not vary across the data set. What the observations resolve is the difference between $\psi$ and the limiting magnitude of each observation, and averaging the transformed magnitudes before converting them back estimates that difference at the reference point. One qualification belongs here. When the limiting magnitudes differ, the mean of the transformed magnitudes is not exactly the value that belongs to their mean, because the back-transformation is curved rather than straight. The resulting shift is small, and it shrinks the more similar the limiting magnitudes are. If all observations share the same limiting magnitude, no averaging over differing limiting magnitudes is involved at all, and the reference point is exact. By contrast, if $\psi$ is significantly larger than the limiting magnitude, estimation becomes problematic: $\psi$ effectively tends to infinity, and the observable magnitude distribution approaches the geometric model of visual meteor magnitudes with a population index of $r \approx 2.5$. In practice, however, it is preferable to use a confidence interval estimate. For example, one can estimate that $\psi$ is, with 10 percent probability, not smaller than: ```{r, echo=TRUE} print(vmideal_vst_to_psi(qnorm(0.90, tm_mean$val, tm_mean$sd), lim_magn_mean)) ``` How accurate the result is depends on two separate things. Over the range the transformation covers — up to about three magnitudes above the limiting magnitude — it reproduces $\psi$ to within a few hundredths; this is the error of the back-transformation, and it stays well below the sampling error even for very large data sets, so it is not what limits an estimate. What does limit it is that sampling error, and here the transformation is somewhat less precise than the maximum likelihood estimate of the previous section. What it offers in exchange is speed. The transformation is applied once per meteor and everything afterwards is an ordinary mean, whereas maximizing the likelihood iterates over the distribution. For a single estimate this hardly matters, as both are fast enough in absolute terms. It becomes decisive as soon as the estimate is repeated: the bootstrap above draws fifty thousand replicates and still runs in a fraction of a second. The same applies when $\psi$ is estimated separately for many time intervals or showers. It is tempting to go one step further. Since the transformed magnitudes are homoscedastic, they look like the ideal input for `lm()`, with the limiting magnitude and any covariate of interest entering as regressors. That step does not work, and it is worth saying why. The expected value of `tm` depends on $\psi$ and the limiting magnitude only through their difference, which is what makes a single reference point meaningful in the first place. That dependence, however, is not linear: it falls off roughly exponentially. A straight line on the `tm` scale is therefore a curve on the $\psi$ scale, and a $\psi$ profile that is in truth linear in a covariate comes back out of such a fit bent. The estimate is biased, not merely imprecise. What makes this worth a warning is that the usual diagnostics do not reveal it. Precisely because the transformation makes the residuals homoscedastic, such a fit looks unremarkable: the residual variance is where it should be, and tests for curvature rarely respond. The discrepancy sits in the expected value rather than in the dispersion, so a clean residual plot is no evidence that the model holds here. Nor does a change of scale help — taking logarithms, for instance, improves the shape but not the conclusion, because the influence of the limiting magnitude then depends on the very $\psi$ that is being estimated. The geometric model behaves differently in this respect, which is instructive. Its transformation is built on the ratio of the perception probabilities of two adjacent magnitude classes, whose mean is $1/r$ and, unlike the magnitudes themselves, no longer depends on the limiting magnitude at all. Accordingly `vmgeom_vst_to_r()` takes no limiting magnitude: there is no reference point to fix, and none to carry through a model. That is what makes the transformed magnitudes usable with `lm()` there, and it is precisely what the ideal distribution does not offer, since the limiting magnitude cannot be eliminated from $\psi$ but only measured against it. The difference also shows in the form of the two back-transformations. For the geometric model it is a power relation, hence an exact straight line on the logarithmic scale, so that effects combine additively there. Here it is a calibrated polynomial in $\log (\texttt{tm})$ with no such form, and no change of scale produces one. Of the three properties by which transformations of this kind are classically judged — constant variance, symmetry and additivity — this one delivers the first and not the last. Beyond estimating $\psi$, the common scale is useful in its own right. Since transformed magnitudes carry the same variance regardless of the limiting magnitude they were observed at, subsets of the data can be compared with each other directly — the mean of `tm` per observer or per night, for instance — which is not possible on the raw magnitude scale, where any such comparison is confounded by the limiting magnitude. For a $\psi$ that varies with a covariate, neither this transformation nor the maximum likelihood estimate of the previous section is sufficient. That is the subject of the next section. ## Generalized Linear Model Both methods so far assume that $\psi$ is a single constant. A generalized linear model lifts that restriction: it estimates $\psi$ as a function of covariates, and the limiting magnitude is carried along per observation instead of being reduced to a common reference point. `vmideal_glm()` fits the ideal distribution through `glm()`. As for the geometric model, the response is passed as a two-column matrix holding the meteor magnitude and the corresponding limiting magnitude, so that the limiting magnitude 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 <- vmideal_glm( cbind(magn, lim_magn) ~ 1, data = magnitudes_obs, weights = magnitudes_obs$Freq ) ``` The link is the identity, since $\psi$ is unrestricted on the real line. `predict()` therefore returns $\psi$ itself. ```{r, echo=TRUE} # the model has no covariates, so any single row predicts the global psi newdata <- data.frame(row.names = "") print(c( "ML" = result_ml$par, "GLM" = as.numeric(predict(result_glm, newdata)) )) ``` Unlike `vmgeom_glm()`, the two estimates do not agree exactly, and this is a property of the model rather than of the implementation. A `glm()` fit works entirely through the expected magnitude $\mu = \mathrm{E}[M]$, which the model predicts for a given limiting magnitude and parameter. It compares each observed magnitude $m$ with that expectation and adjusts the parameter by the difference $m - \mu$. This reproduces the maximum likelihood estimate exactly whenever the deviation $m - \mu$ carries all the information the magnitude has about the parameter — in other words, whenever the mean is a *sufficient statistic*. For the geometric model that is the case: $\log f$ is linear in the magnitude, so a meteor twice as far from the expectation contributes twice as much. The density of the ideal distribution contains the term $r^\psi + r^m$, which makes $\log f$ non-linear in $m$. How much a faint meteor really tells us about $\psi$ therefore levels off, while the difference $m - \mu$ that `glm()` uses in its place keeps growing. The two stop being proportional, and the mean is not sufficient for $\psi$. `vmideal_glm()` is therefore a quasi-likelihood estimator. It is consistent, but not fully efficient, and how much it loses depends on how far $\psi$ lies above the limiting magnitude. While $\psi$ stays below it, the variance is within about ten percent of the Cramér-Rao bound; once $\psi$ exceeds it, the factor grows to roughly $1.4$ to $2.2$, the larger values belonging to the smaller limiting magnitudes. The difference of more than a magnitude between the two estimates above is a case in point. This excerpt contains only `r sum(magnitudes$Freq)` meteors at limiting magnitudes between `r min(magnitudes$lim_magn)` and `r max(magnitudes$lim_magn)`, and the estimate already sits above the mean limiting magnitude, where the mean of the distribution begins to flatten. The two estimators are then far less tightly coupled than they are for a few thousand meteors at a moderate $\psi$, where they typically differ by a few hundredths. The standard error reflects this honestly: ```{r, echo=TRUE} print(predict(result_glm, newdata, se.fit = TRUE)) ``` That is also the limit worth knowing. What the data can resolve is governed by the distance between $\psi$ and the limiting magnitude, not by $\psi$ itself. Up to about two to three magnitudes above the limiting magnitude $\psi$ is estimated reliably; beyond that the mean of the distribution flattens, and once it has reached the value it converges to, no data can tell a larger $\psi$ from a smaller one. This never fails silently. The standard error grows as the mean flattens, and where $\psi$ is no longer identified at all, `predict()` returns `Inf`: ```{r, echo=TRUE} # magnitudes fainter than the ideal distribution can produce faint <- data.frame(magn = c(rep(3L, 40), rep(4L, 40), rep(5L, 20)), lim_magn = 5.5) result_faint <- vmideal_glm(cbind(magn, lim_magn) ~ 1, data = faint) psi_faint <- predict(result_faint, faint[1, ]) print(psi_faint) ``` `Inf` is a result rather than an error: it says the magnitudes are as faint as the model can make them, so $\psi$ is bounded below by the fit rather than determined by it. Such data are either a small-sample accident or a sign that the ideal distribution is not the right model for them. The linear predictor stays finite in both cases, so `summary()` and `anova()` remain usable. Being a result, it can be used as one. `dvmideal()` and its companions accept an infinite $\psi$ and evaluate the distribution the ideal one converges to, so the prediction needs no special handling before it is passed on: ```{r, echo=TRUE} print(dvmideal(seq(5L, 1L), 5.5, psi_faint)) ``` Those are the probabilities of the geometric model with a population index of $r = 10^{0.4}$, which is what an unbounded $\psi$ means: ```{r, echo=TRUE} print(identical( dvmideal(seq(5L, 1L), 5.5, psi_faint), dvmgeom(seq(5L, 1L), 5.5, 10^0.4) )) ``` Just short of that point the mean moves so little with $\psi$ that the iteration can run out of steps. `vmideal_glm()` then warns instead of returning the value it happened to stop at; the estimate is unreliable there, and its standard error does not show it. For a single global $\psi$ the two methods of the previous sections remain preferable — maximum likelihood for accuracy, the transformation for speed. The benefit of the linear model 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 $\psi$ varies with solar longitude would be written as ```r vmideal_glm(cbind(magn, lim_magn) ~ sl, data = magnitudes, weights = Freq) ``` and a smooth activity-dependent profile as ```r vmideal_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()` is not restricted to comparing models of the same kind. The geometric model of visual meteor magnitudes can be fitted to the same data with `vmgeom_glm()` and compared directly: ```{r, echo=TRUE} result_geom <- vmgeom_glm( cbind(magn, lim_magn) ~ 1, data = magnitudes_obs, weights = magnitudes_obs$Freq ) print(c("ideal" = AIC(result_glm), "geometric" = AIC(result_geom))) ``` Both models have one parameter here, so the comparison reduces to their likelihoods. The lower value belongs to the ideal distribution, but the two lie close enough together that this excerpt of `r sum(magnitudes$Freq)` meteors does not decide between them; a difference of a few units is weak evidence. The comparison is only meaningful because both fits use the same observations and the same response — a model fitted to a subset, or with the zero-frequency rows retained, is not comparable in this way. ## 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. The maximum likelihood estimate of the first section is used for this. What the test is meant to reveal is whether the model itself describes the data, and any inefficiency of the estimator would enter the result alongside it. Since the maximum likelihood estimate is the one that extracts all the information the magnitudes carry about $\psi$, a departure found here can be attributed to the model rather than to the way the parameter was obtained. The correction of the degrees of freedom applied below rests on the same assumption. Using the estimated parameter, we retrieve the relative frequencies `p` for each observation and add them to the data frame `magnitudes`: ```{r, echo=TRUE, results='asis'} psi_mean <- result_ml$par magnitudes$p <- with(magnitudes, dvmideal(m = magn, lm = lim_magn, psi_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 `pvmideal()` to determine their probability. ```{r, echo=TRUE, results='asis'} idx <- magnitudes$magn == magn_min magnitudes$p[idx] <- with( magnitudes[idx, ], pvmideal(m = magn + 1L, lm = lim_magn, psi_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 the ideal meteor magnitude distribution and that the perception probabilities are correct (i.e., error-free or precisely known), these 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)) ```