Подходящее распределение к заданным значениям частоты в R

Вот более удачная попытка, например, до использования optim для поиска наилучшего значения, ограниченного набор значений в поле (определяется векторами lower и upper в вызове optim). Обратите внимание, что он масштабирует x и y как часть оптимизации в дополнение к параметру формы распределения Вейбулла, поэтому у нас есть 3 параметра для оптимизации.

К сожалению, при использовании всех точек он почти всегда находит что-то на краях ограничивающей рамки, что указывает мне на то, что, возможно, Вейбулл не подходит для всех данных. Проблема в двух точках - они слишком большие. Вы видите попытку подгонки ко всем данным на первом графике.

Если я отброшу эти первые два пункта и просто подгоню остальные, мы получим гораздо лучшее соответствие. Вы видите это на втором графике. Я думаю, что это хорошее соответствие, в любом случае это локальный минимум внутри ограничивающей коробки.

library(optimx)
sample <- c(60953,7787,3056,2359,1759,1819,1189,1077,1080,985,622,648,518,
            611,1037,727,489,432,371,1125,69,595,624)
t.sample <- 0:22

s.fit <- sample[3:23]
t.fit <- t.sample[3:23]

wx <- function(param) { 
  res <- param[2]*dweibull(t.fit*param[3],shape=param[1])
  return(res)
} 
minwx <- function(param){
  v <- s.fit-wx(param)
  sqrt(sum(v*v))
}

p0 <- c(1,200,1/20)
paramopt <- optim(p0,minwx,gr=NULL,lower=c(0.1,100,0.01),upper=c(1.1,5000,1))

popt <- paramopt$par
popt
rms <- paramopt$value
tit <- sprintf("Weibull - Shape:%.3f xscale:%.1f  yscale:%.5f rms:%.1f",popt[1],popt[2],popt[3],rms)

plot(t.sample[2:23], sample[2:23], type = "p",col="darkred")
lines(t.fit, wx(popt),col="blue")
title(main=tit)

Вы можете напрямую рассчитать параметры максимальной вероятности, как описано здесь.

# Defining the error of the implicit function
k.diff <- function(k, vec){
  x2 <- seq(length(vec))
  abs(k^-1+weighted.mean(log(x2), w = sample)-weighted.mean(log(x2), 
                                                            w = x2^k*sample))
}

# Setting the error to "quite zero", fulfilling the equation
k <- optimize(k.diff, vec=sample, interval=c(0.1,5), tol=10^-7)$min

# Calculate lambda, given k
l <- weighted.mean(seq(length(sample))^k, w = sample)

# Plot
plot(density(rep(seq(length(sample)),sample)))
x <- 1:25
lines(x, dweibull(x, shape=k, scale= l))

Предполагая, что данные получены из распределения Вейбулла, вы можете получить оценку параметра формы и масштаба следующим образом:

sample <- c(7787,3056,2359,1759,1819,1189,1077,1080,985,622,648,518,
        611,1037,727,489,432,371,1125,69,595,624)
 f<-fitdistr(sample, 'weibull')
 f

Если вы не уверены, распространяется ли Weibull, я бы рекомендовал использовать ks.test. Это проверяет, являются ли ваши данные гипотетическим распределением. Зная ваши знания о характере данных, вы можете протестировать несколько выбранных дистрибутивов и посмотреть, какой из них работает лучше всего.

Для вашего примера это будет выглядеть так:

 ks = ks.test(sample, "pweibull", shape=f$estimate[1], scale=f$estimate[2])
 ks

Значение p незначительно, поэтому вы не отвергаете гипотезу о том, что данные получены из распределения Вейбулла.

Обновление: гистограммы Вейбулла или экспоненты хорошо соответствуют вашим данным. Я думаю, что экспоненциальное распределение дает вам лучшее соответствие. Распределение Парето — еще один вариант.

f<-fitdistr(sample, 'weibull')
z<-rweibull(10000, shape= f$estimate[1],scale= f$estimate[2])
hist(z)

f<-fitdistr(sample, 'exponential')
z = rexp(10000, f$estimate[1]) 
hist(z)