使用mmedist或fitdist（with mme）误差估计Frechet分布的参数

问题描述

我在R语言方面相对较新，如果您可以看一下以下代码，将不胜感激。我正在尝试使用mmedist估算Frechet分布（或反韦布尔）的形状参数（我也尝试了调用mmedist的fitdist），但似乎出现以下错误：

I'm relatively new in R and I would appreciated if you could take a look at the following code. I'm trying to estimate the shape parameter of the Frechet distribution (or inverse weibull) using mmedist (I tried also the fitdist that calls for mmedist) but it seems that I get the following error :

mmedist中的错误（数据，distname，start = start，fix.arg = fix.arg，...）：
必须定义经验矩函数。 / code>

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

请注意，我多次尝试更改代码，以查看我的错误是否在语法上但我总是会遇到相同的错误。

Please note that I tried many times en-changing the code in order to see if my mistake was in syntax but I always get the same error.

我知道文档中的范例。我也尝试过，但是没有运气。请注意，为了使该方法起作用，力矩的顺序必须小于shape参数（即shape）。

I'm aware of the paradigm in the documentation. I've tried that as well but with no luck. Please note that in order for the method to work the order of the moment must be smaller than the shape parameter (i.e. shape).

示例如下：

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)

谢谢

推荐答案

您将需要对 mmedist -我建议您复制代码，并创建自己的函数。

You will need to make non-trivial changes to the source of mmedist -- I recommend that you copy out the code, and make your own function foo_mmedist.

1. 您需要进行的第一个更改是在 mmedist 的第94行：

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

该行检查 memp是否存在，而不是是否存在您实际传递的参数作为函数存在。

That line checks whether "memp" is a function that exists, as opposed to whether the argument that you have actually passed exists as a function.

if（！exists（as.character（expression（memp）），mode = function））

第二点，正如我已经指出的，与 optim 例程实际调用的事实有关 funobj 调用 DIFF2 ，后者调用（参见第112行）用户提供的 memp 函数，在您的情况下为 minvweibull ，带有两个参数- obs ，其解析为数据和订单，但是由于 minvweibull 不会将数据作为第一个参数，它将失败。

这是预期的，因为帮助页面会告诉您：

The second, as I have already noted, relates to the fact that the optim routine actually calls funobj which calls DIFF2, which calls (see line 112) the user-supplied memp function, minvweibull in your case with two arguments -- obs, which resolves to data and order, but since minvweibull does not take data as the first argument, this fails.
This is expected, as the help page tells you:

如何解决此问题？从 moments 包传递函数 moment 。这是完整的代码（假设您已进行上述更改，并创建了一个名为 foo_mmedist 的新函数）：

How can you fix this? Pass the function moment from the moments package. Here is complete code (assuming that you have made the change above, and created a new function called 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)

您可以检查优化是否按预期进行：

You can check that optimization has occurred as expected:

> para_lm$estimate
   shape    scale
2.490816 1.004128

但是请注意，实际上简化为一种过分确定的矩量法的粗略方法，并且不确定这在理论上是否合适。

Note however, that this actually reduces to a crude way of doing overdetermined method of moments, and am not sure that this is theoretically appropriate.

这篇关于使用mmedist或fitdist（with mme）误差估计Frechet分布的参数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！