mideal {vismeteor} | R Documentation |
Density, distribution function, quantile function and random generation of ideal distributed meteor magnitudes.
dmideal(m, psi = 0, log = FALSE)
pmideal(m, psi = 0, lower.tail = TRUE, log = FALSE)
qmideal(p, psi = 0, lower.tail = TRUE)
rmideal(n, psi = 0)
m |
numeric; meteor magnitude. |
psi |
numeric; the location parameter of a probability distribution. It is the only parameter of the distribution. |
log |
logical; if |
lower.tail |
logical; if |
p |
numeric; probability. |
n |
numeric; count of meteor magnitudes. |
The density of an ideal magnitude distribution is
{\displaystyle \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
is the meteor magnitude, r = 10^{0.4} \approx 2.51189 \dots
is a constant and
\psi
is the only parameter of this magnitude distribution.
dmideal
gives the density, pmideal
gives the distribution function,
qmideal
gives the quantile function and rmideal
generates random deviates.
The length of the result is determined by n
for rmideal
, and is the maximum
of the lengths of the numerical vector arguments for the other functions.
qmideal
can return NaN
value with a warning.
Richter, J. (2018) About the mass and magnitude distributions of meteor showers. WGN, Journal of the International Meteor Organization, vol. 46, no. 1, p. 34-38
old_par <- par(mfrow = c(2,2))
psi <- 5.0
plot(
function(m) dmideal(m, psi, log = FALSE),
-5, 10,
main = paste0('density of ideal meteor magnitude\ndistribution (psi = ', psi, ')'),
col = "blue",
xlab = 'm',
ylab = 'dp/dm'
)
abline(v=psi, col="red")
plot(
function(m) dmideal(m, psi, log = TRUE),
-5, 10,
main = paste0('density of ideal meteor magnitude\ndistribution (psi = ', psi, ')'),
col = "blue",
xlab = 'm',
ylab = 'log( dp/dm )'
)
abline(v=psi, col="red")
plot(
function(m) pmideal(m, psi),
-5, 10,
main = paste0('probability of ideal meteor magnitude\ndistribution (psi = ', psi, ')'),
col = "blue",
xlab = 'm',
ylab = 'p'
)
abline(v=psi, col="red")
plot(
function(p) qmideal(p, psi),
0.01, 0.99,
main = paste('quantile of ideal meteor magnitude\n distribution (psi = ', psi, ')'),
col = "blue",
xlab = 'p',
ylab = 'm'
)
abline(h=psi, col="red")
# generate random meteor magnitudes
m <- rmideal(1000, psi)
# log likelihood function
llr <- function(psi) {
-sum(dmideal(m, psi, log=TRUE))
}
# maximum likelihood estimation (MLE) of psi
est <- optim(2, llr, method='Brent', lower=0, upper=8, hessian=TRUE)
# estimations
est$par # mean of psi
sqrt(1/est$hessian[1][1]) # standard deviation of psi
par(old_par)