2013-03-08 2 views
1

Я использовал функцию fitdistr из пакета R MASS для настройки функции плотности вероятности параметров Weibull 2 ​​(pdf).Функция кумулятивного распределения Вейбула, начиная с команды «fitdistr»

Это мой код:

require(MASS) 

h = c(31.194, 31.424, 31.253, 25.349, 24.535, 25.562, 29.486, 25.680, 26.079, 30.556,  30.552, 30.412, 29.344, 26.072, 28.777, 30.204, 29.677, 29.853, 29.718, 27.860, 28.919, 30.226, 25.937, 30.594, 30.614, 29.106, 15.208, 30.993, 32.075, 31.097, 32.073, 29.600, 29.031, 31.033, 30.412, 30.839, 31.121, 24.802, 29.181, 30.136, 25.464, 28.302, 26.018, 26.263, 25.603, 30.857, 25.693, 31.504, 30.378, 31.403, 28.684, 30.655, 5.933, 31.099, 29.417, 29.444, 19.785, 29.416, 5.682, 28.707, 28.450, 28.961, 26.694, 26.625, 30.568, 28.910, 25.170, 25.816, 25.820) 

weib = fitdistr(na.omit(h),densfun=dweibull,start=list(scale=1,shape=5)) 

hist(h, prob=TRUE, main = "", xlab = "x", ylab = "y", xlim = c(0,40), breaks = seq(0,40,5)) 
curve(dweibull(x, scale=weib$estimate[1], shape=weib$estimate[2]),from=0, to=40, add=TRUE) 

Теперь я хотел бы создать интегральную функцию распределения Вейбулла (КОР) и построить его в виде графика:

enter image description here, где х> 0, б = шкала, a = форма

Я попытался применить параметры шкалы и формы для h, используя приведенную выше формулу, но это было не так.

ответ

3

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

cdweibull <- function(x, shape, scale, log = FALSE){ 
    dd <- dweibull(x, shape= shape, scale = scale, log = log) 
    dd <- cumsum(dd) * c(0, diff(x)) 
    return(dd) 
} 

Обсуждение выше о масштабе различий несмотря на это, вы можете построить его над графиком же, как dweibull:

require(MASS) 

h = c(31.194, 31.424, 31.253, 25.349, 24.535, 25.562, 29.486, 25.680, 
     26.079, 30.556, 30.552, 30.412, 29.344, 26.072, 28.777, 30.204, 
     29.677, 29.853, 29.718, 27.860, 28.919, 30.226, 25.937, 30.594, 
     30.614, 29.106, 15.208, 30.993, 32.075, 31.097, 32.073, 29.600, 
     29.031, 31.033, 30.412, 30.839, 31.121, 24.802, 29.181, 30.136, 
     25.464, 28.302, 26.018, 26.263, 25.603, 30.857, 25.693, 31.504, 
     30.378, 31.403, 28.684, 30.655, 5.933, 31.099, 29.417, 29.444, 
     19.785, 29.416, 5.682, 28.707, 28.450, 28.961, 26.694, 26.625, 
     30.568, 28.910, 25.170, 25.816, 25.820) 

weib = fitdistr(na.omit(h),densfun=dweibull,start=list(scale=1,shape=5)) 

hist(h, prob=TRUE, main = "", xlab = "x", 
    ylab = "y", xlim = c(0,40), breaks = seq(0,40,5), ylim = c(0,1)) 

curve(cdweibull(x, scale=weib$estimate[1], shape=weib$estimate[2]), 
    from=0, to=40, add=TRUE) 

histogram with weibull cumulative density overlay

+0

Он может работать для точек выборки с равномерным интервалом больше одного, но график, по-видимому, менее симпатичный. – Noah

+0

Это здорово! Просто одно сомнение. Настройка, о которой вы говорили, находится в аргументе 'lag' внутри' diff (x) '? Например: значение по умолчанию для 'lag' равно 1, но если бы у моего образца было большее расстояние, я бы установил lag = 2,4,10 и т. Д.. Разве это так? –

+0

Andre, я так не думаю, потому что 'cumsum' работает непосредственно над результатами, которые возвращаются' dweibull'. Коррекция «отставания», в некотором смысле, является _what возвращается 'diff'. Вы можете попробовать: после запуска команды 'hist' выше, запустите' lines (seq (0, 40, by = 4), cdweibull (seq (0, 40, by = 4), scale = weib $ оценка [1] , shape = weib $ оценка [2])). Вы увидите, что строка добавляется, но она более грубая, потому что разрешение намного ниже. – Noah

3

Это работает для моих данных, но ваше может отличаться.

>h=rweibull3(1000,2,2,2) 

>#this gives some warnings...that I ignore. 
>weib = fitdistr(h,densfun=dweibull3,start=list(scale=1,shape=5,thres=0.5)) 

There were 19 warnings (use warnings() to see them)  

Означает, что начальные значения влияют на то, как происходит приход. Поэтому, если начальные значения близки к истинным значениям, вы получите меньше предупреждений.

>curve(dweibull3( x, 
      scale=weib$estimate[1], 
      shape=weib$estimate[2], 
      thres=weib$estimate[3]), 
      add=TRUE) 

enter image description here

+0

ли Вейбулла плотность определяется для диапазона значений, которые вы хотите, чтобы соответствовать? Он не работает для значений, равных нулю или меньше. – Seth

Смежные вопросы