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

У меня меняется значение частоты с единицами времени (x), как показано на рисунке ниже. После некоторой нормализации эти значения можно рассматривать как точки данных функции плотности для некоторого распределения.

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

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

plot(1:length(sample), sample, type = "l")
points(1:length(sample), sample)

enter image description here

Обновление. Чтобы не допустить неправильного понимания, я хотел бы добавить немного больше объяснений. Говоря, что меняют значения частоты, меняющиеся с помощью единиц времени (x), я имею в виду, что у меня есть данные, которые говорят, что у меня есть:

  • 7787 реализаций значения 1
  • 3056 реализации значения 2
  • 2359 реализация значений 3... и т.д.

Какой-то путь к моей цели (как мне кажется, некорректный) - создать набор этих реализаций:

# Loop to simulate values 
set.values <- c()
for(i in 1:length(sample)){
  set.values <<- c(set.values, rep(i, times = sample[i]))
}

hist(set.values)
lines(1:length(sample), sample)
points(1:length(sample), sample)

enter image description here

и используйте fitdistr на set.values:

f2 <- fitdistr(set.values, 'weibull')
f2

Почему я думаю, что это неправильный путь и почему я ищу лучшее решение в R?

  • в представленном выше подходе к распределению, предполагается, что set.values является полным набором моих реализаций из распределения T

  • в моем первоначальном вопросе я знаю точки из первой части кривой плотности - я не знаю его хвоста, и я хочу оценить хвост (и всю функцию плотности)

Ответы

Ответ 1

First try with all points

Second try with first point dropped Вот лучшая попытка, как прежде, чем использовать optim, чтобы найти наилучшее значение, ограниченное набором значений в поле (определяемом векторами lower и upper в вызове optim). Обратите внимание, что он масштабирует x и y как часть оптимизации в дополнение к параметру формы распределения Weibull, поэтому у нас есть 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)

Ответ 2

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

# 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))

Ответ 3

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

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)