Оценка моментов от функции распределения - PullRequest
Новая
       90

Оценка моментов от функции распределения

0 голосов
/ 30 июня 2018

У меня есть ненормальная функция распределения, которую я хочу вычислить по ее моментам (среднее значение, дисперсия, асимметрия и эксцесс). Я знаю пакеты e1071 и moments, которые вычисляют моменты для вектора дискретных значений. Есть ли пакет, который оценивает моменты для непрерывной функции распределения? В качестве примера предположим, что у меня есть этот дистрибутив: 

tmp <- runif(100,-10,10)
PDFdata <- density(tmp , kernel = "gaussian")
PDF <- splinefun(PDFdata$x , PDFdata$y)

Теперь я хочу вычислить:

enter image description here

Ответы [ 2 ]

0 голосов
/ 30 июня 2018

Если вы оцениваете свою плотность по данным, вам лучше использовать эмпирические моменты из данных для оценки моментов распределения. Если вы просто использовали это как пример функции, то вы могли бы использовать функцию integrate из пакета stats. Например, 

set.seed(123)
tmp <- runif(100,-10,10)
PDFdata <- density(tmp , kernel = "gaussian")
PDF <- splinefun(PDFdata$x , PDFdata$y)

mean(tmp)
[1] -0.02882012
integrate(function(x) x*PDF(x), -Inf, Inf)
Error in integrate(function(x) x * PDF(x), -Inf, Inf) : 
  the integral is probably divergent

и действительно, для этого набора данных функция PDF не является плотностью: 

plot(PDF, from = -100, to = 100)

enter image description here

Таким образом, мы заставляем его обнуляться вне диапазона оценки плотности: 

PDF2 <- function(x) ifelse(x < min(PDFdata$x) | x > max(PDFdata$x), 0, 
                           PDF(x))

integrate(function(x) x*PDF2(x), -Inf, Inf)
-0.02900585 with absolute error < 8.2e-05
0 голосов
/ 30 июня 2018

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

Во-первых, я сделаю результаты воспроизводимыми, установив начальное значение ГСЧ. 

set.seed(6126)    # Make the results reproducible

tmp <- runif(100,-10,10)
PDFdata <- density(tmp , kernel = "gaussian")
PDF <- splinefun(PDFdata$x , PDFdata$y)

int <- function(x, c = 0, g = PDF, n = 1) (x - c)^n * g(x)

integrate(int, lower = -15, upper = 15, n = 1)
#-0.3971095 with absolute error < 7.8e-06

integrate(int, lower = -15, upper = 15, n = 2)
#35.76295 with absolute error < 0.0012

plot(PDFdata)
curve(PDF, -15, 15, add = TRUE, col = "blue", lty = "dashed")

density plot

...