Насколько я понимаю, формула представляет собой вероятность обнаружения молекулы на определенном расстоянии D * после n шагов, с учетом константы диффузии D . Поскольку D * немного сбивает с толку, потому что его можно перепутать с D , он был заменен в функции на x . Обратите внимание, что в формуле (n - 1)! - это просто гамма-функция gamma(n). Следовательно, функция в R будет записана как:
f <- function(x, n, D) (n/D)^n * x^(n-1) * exp(-n * x/D) / gamma(n)
Проблема в том, что, хотя данные содержат значения, которые мы можем подставить для x и n (df $ D и df $ n), у нас нет константы диффузии, и на самом деле это, по-видимому, то, что вы пытаетесь здесь найти.
Однако мы можем использовать указанную выше функцию, чтобы получить сумму отрицательных логарифмических вероятностей наличия эти точки данных для любого заданного значения D
f2 <- function(D) sum(-log(f(df$D, df$n, D)/(1 - f(df$D, df$n, D))))
Чтобы упростить задачу, мы векторизуем его:
f3 <- function(D) sapply(D, f2)
Это означает, что если мы построим различные значений для D , мы можем найти оценку максимального правдоподобия для D :
plot(0:100, f3(0:100))
Из этого следует, что наше оптимальное значение D находится где-то между 4 и 9:
plot(seq(4, 9, 0.1), f3(seq(4, 9, 0.1)))
и между 7 и 8:
plot(seq(7, 8, 0.01), f3(seq(7, 8, 0.01)))
Мы можем найти точную точку следующим образом:
optimize(f3, range = c(7, 8))
#> $minimum
#> [1] 7.442408
#> $objective
#> [1] 5984.383
Это означает, что наш Лучший предположение для значения D составляет 7,442408 согласно вашим данным. Таким образом, мы можем вычислить плотность вероятности обнаружения молекулы на определенном расстоянии для всех различных значений n :
n <- c(4:10, 20, 100)
x <- seq(0, 30, 0.1)
dens <- do.call(c, lapply(n, function(i) f(x, i, 7.442408)))
sim <- data.frame(x = rep(x, length(n)), y = dens, n = factor(rep(n, each = 301)))
ggplot(sim, aes(x, y, colour = n)) + geom_line() + theme_minimal()
Итак, мы видим, что результатом является семейство гамма-распределений, которые меняются в зависимости от количества шагов.
Чтобы убедить себя, что мы на правильном пути, мы может построить эмпирическую функцию плотности наших точек данных при n == 4 и наложить нашу предсказанную функцию плотности:
plot(density(df$D[df$n == 4]), ylim= c(0, 0.2))
lines(sim$x[sim$n == 4], sim$y[sim$n == 4], lty = 2, col = 2)
Это выглядит довольно близко .