Я думаю, что эта проблема в конечном итоге является проблемой оптимизации, когда обрабатывается одна строка данных за раз. Однако, поскольку вы хотите масштабировать его, вот пример для нахождения параметров ядра распределения.
Этот процесс не оптимизация: он расширяет определенный диапазон возможных параметров k (форма) и находит комбинацию формы / масштаба (с учетом вашего среднего значения), которая наиболее близко соответствует вашим верхним и нижним квантилям. Вы контролируете степень детализации k, которая настолько хороша, насколько вы собираетесь получить допуск (что было бы подходящим для оптимизации).
Таким образом, этот процесс будет несовершенным . Я надеюсь, что это даст вам достаточно быстрый процесс для достаточно хорошей точности.
Я собираюсь сначала продемонстрировать процесс, который обрабатывает одну строку за раз, как guesser1. Затем я разверну его, чтобы выполнить ту же операцию с произвольным числом mean, lower и upper.
Данные с известными ответами
Но сначала я хочу чтобы сгенерировать мои собственные образцы, чтобы мы знали «правду» для этого угадателя.
library(dplyr)
set.seed(42)
n <- 4
randks <- tibble(
k = sample(1:10, size = n, replace = TRUE),
scale = sample(seq(0.5, 2.5, by = 0.5), size = n, replace = TRUE)
) %>%
mutate(
samp = map2(k, scale, ~ rgamma(1000, shape = .x, scale = .y)),
theor_mu = k*scale,
mu = map_dbl(samp, ~ mean(.x)),
lwr = map_dbl(samp, ~ quantile(.x, 0.025)),
upr = map_dbl(samp, ~ quantile(.x, 0.975))
) %>%
select(-samp)
randks
# # A tibble: 4 x 6
# k scale theor_mu mu lwr upr
# <int> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 10 2 20 19.9 9.47 33.7
# 2 10 1.5 15 15.1 7.36 25.0
# 3 3 2 6 5.85 1.08 14.5
# 4 9 0.5 4.5 4.51 1.99 7.72
Guesser1
По одной строке за раз:
guesser1 <- function(mu, lwr, upr, k.max = 10, k.by = 0.01) {
stopifnot(length(mu) == 1, length(lwr) == 1, length(upr) == 1)
ks <- seq(0, k.max, by = k.by)[-1]
L <- sapply(ks, function(k) qgamma(0.025, shape = k, scale = mu / k))
U <- sapply(ks, function(k) qgamma(0.975, shape = k, scale = mu / k))
dists <- sqrt((L-lwr)^2 + (U-upr)^2)
ind <- which.min(dists)
data.frame(
k = ks[ind],
scale = mu / ks[ind],
dist = min(dists),
lwr = L[ind],
upr = U[ind]
)
}
В действии :
out1 <- do.call(rbind, Map(guesser1, randks$mu, randks$lwr, randks$upr))
cbind(subset(randks, select = -theor_mu), out1)
# k scale mu lwr upr k scale dist lwr upr
# 1 10 2.0 19.88 9.47 33.67 10.00 1.988 0.304 9.53 33.97
# 2 10 1.5 15.06 7.36 25.02 10.00 1.506 0.727 7.22 25.73
# 3 3 2.0 5.85 1.08 14.50 2.76 2.120 0.020 1.10 14.50
# 4 9 0.5 4.51 1.99 7.72 9.55 0.472 0.142 2.12 7.79
### \____ randks __________/ \____ guessed ____________/
Определенно есть некоторые различия, подчеркивающие мое первоначальное утверждение, что это несовершенно.
Гадалки
Все строки сразу. Это немного больше работы в функции, поскольку она имеет дело с матрицами, а не только с векторами. Не проблема, я просто хотел доказать это по одному, прежде чем отправиться на поиски удовольствия.
guessers <- function(mu, lwr, upr, k.max = 10, k.by = 0.01, include.size = FALSE) {
stopifnot(length(mu) == length(lwr), length(mu) == length(upr))
# count <- length(mu)
ks <- seq(0, k.max, by = k.by)[-1]
# 'ks' dims: [mu]
L <- outer(mu, ks, function(m, k) qgamma(0.025, shape = k, scale = m / k))
U <- outer(mu, ks, function(m, k) qgamma(0.975, shape = k, scale = m / k))
# 'L' & 'U' dims: [mu, ks]
dists <- sqrt((L - lwr)^2 + (U - upr)^2)
inds <- apply(dists, 1, which.min)
mindists <- apply(dists, 1, min)
i <- seq_along(mu)
out <- data.frame(
k = ks[inds],
scale = mu / ks[inds],
dist = mindists,
lwr = L[cbind(i, inds)],
upr = U[cbind(i, inds)]
)
size <- if (include.size) {
message("guessers memory: ",
object.size(list(ks, L, U, dists, inds, mindists, i, out)))
}
out
}
В действии:
outs <- guessers(randks$mu, randks$lwr, randks$upr, include.size = TRUE)
# guessers memory: 106400
cbind(subset(randks, select = -theor_mu), outs)
# k scale mu lwr upr k scale dist lwr upr
# 1 10 2.0 19.88 9.47 33.67 10.00 1.988 0.304 9.53 33.97
# 2 10 1.5 15.06 7.36 25.02 10.00 1.506 0.727 7.22 25.73
# 3 3 2.0 5.85 1.08 14.50 2.76 2.120 0.020 1.10 14.50
# 4 9 0.5 4.51 1.99 7.72 9.55 0.472 0.142 2.12 7.79
### \____ randks __________/ \____ guessed (same) _____/
(Я включил сообщение памяти там, чтобы отследить, насколько это может масштабироваться. Сейчас это неплохо, и этот аргумент определенно не должен использоваться в производстве. К вашему сведению.)
Сравнение
Используя microbenchmark, мы повторяем каждую операцию несколько раз и сравните время их выполнения.
microbenchmark::microbenchmark(
g1 = Map(guesser1, randks$mu, randks$lwr, randks$upr),
gs = guessers(randks$mu, randks$lwr, randks$upr)
)
# Unit: milliseconds
# expr min lq mean median uq max neval
# g1 27.3 28.8 33.9 29.7 33.0 131.1 100
# gs 13.3 13.6 14.4 13.8 14.6 20.3 100
Неудивительно, что all-at-once guessers работает немного быстрее. Когда этого не будет? Когда количество строк становится настолько большим, что потребление памяти становится проблемой. Я не знаю, что это такое.
Давайте попробуем то же самое, но увеличим randks с 4 строк до 1000 и повторим тест.
n <- 1000
# randks <- ...
nrow(randks)
# [1] 1000
microbenchmark::microbenchmark(
g1 = Map(guesser1, randks$mu, randks$lwr, randks$upr),
gs = guessers(randks$mu, randks$lwr, randks$upr),
times = 10
)
# Unit: seconds
# expr min lq mean median uq max neval
# g1 8.50 8.99 9.59 9.31 9.64 11.72 10
# gs 3.35 3.44 3.61 3.63 3.77 3.93 10
Так что это определенно быстрее. Среднее время выполнения для 1000 оценок составляет 3,63 секунды, поэтому кажется, что он составляет sh около 300 / se c.
guessers memory: 24066176
(24 MiB) На самом деле, это неплохо для все. Уменьшите k.by, чтобы повысить точность, или увеличьте k.by, чтобы ускорить это.