Оценить параметры распределения Frechet с использованием ошибки mmedist или fitdist (с mme)

Question

Я относительно новичок в R, и я бы оценил, можете ли вы взглянуть на следующий код. Я пытаюсь оценить параметр формы распределения Фреше (или обратный Вейбулла) с помощью mmedist (я попытался также fitdist, который призывает к mmedist), но мне кажется, что я получаю следующее сообщение об ошибке:

Error in mmedist(data, distname, start = start, fix.arg = fix.arg, ...) : the empirical moment function must be defined.

код, который я использую ниже:

require(actuar) 
library(fitdistrplus) 
library(MASS) 
#values 
n=100 
scale = 1 
shape=3 
# simulate a sample 
data_fre = rinvweibull(n, shape, scale) 
memp=minvweibull(c(1,2), shape=3, rate=1, scale=1) 
# estimating the parameters 
para_lm = mmedist(data_fre,"invweibull",start=c(shape=3,scale=1),order=c(1,2),memp = "memp") 

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

Я знаю парадигму в документации. Я тоже пробовал это, но не повезло. Обратите внимание, что для того, чтобы метод работал, порядок момента должен быть меньше, чем параметр формы (т. Е. Форма).

Примером является следующее:

require(actuar) 
#simulate a sample 
x4 <- rpareto(1000, 6, 2) 

#empirical raw moment 
memp <- function(x, order) 
    ifelse(order == 1, mean(x), sum(x^order)/length(x)) 


#fit 
mmedist(x4, "pareto", order=c(1, 2), memp="memp", 
    start=c(shape=10, scale=10), lower=1, upper=Inf) 

Заранее спасибо за любую помощь.

Answer 1

Вам необходимо будет сделать нетривиальные изменения в источнике mmedist. Я рекомендую вам скопировать код и сделать свою собственную функцию foo_mmedist.

1. Первое изменение, которое вы должны сделать это на линии 94 из mmedist:

   if (!exists("memp", mode = "function")) 
   

   Это проверяет линии ли «memp» является функцией, которая существует, в отличие от того аргумента, что у вас есть фактически передается как функция.

   if (!exists(as.character(expression(memp)), mode = "function"))

2. Во-вторых, как уже отмечалось, связано с тем, что optim процедура фактически вызывает funobj который вызывает DIFF2, который вызывает (см линии 112) на предоставленный пользователем memp функцию, minvweibull в ваш случай с двумя аргументами - obs, который разрешается до data и order, но так как minvweibull не принимает данные в качестве первого аргумента, это терпит неудачу.
   Ожидается, что страница справки говорит вам:

   memp Функция реализация эмпирических моментов, сырые или по центру, но должен соответствовать РАСПРУ аргумент. Эта функция должна иметь два аргумента: в качестве первого - числовой вектор данных и как второй порядок момента, возвращаемого функцией.

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

# values 
n = 100 
scale = 1 
shape = 3 

# simulate a sample 
data_fre = rinvweibull(n, shape, scale) 

# estimating the parameters 
para_lm = foo_mmedist(data_fre, "invweibull", 
        start= c(shape=5,scale=2), order=c(1, 2), memp = moment) 

Вы можете проверить, что оптимизация произошло, как и ожидалось:

> para_lm$estimate 
    shape scale 
2.490816 1.004128 

Примечание однако, что это фактически сводится к грубому способу делать переопределенный метод моментов, и я не уверен, что это теоретически подходит.