【问题标题】:Use function mle2() in bbmle package(R) to get parameters for exponential and power law distribution [closed]使用 bbmle 包(R)中的函数 mle2()来获取指数和幂律分布的参数 [关闭]
【发布时间】:2015-12-02 02:01:17
【问题描述】:

这是mydata的一部分:

原始数据很大,所以我上传了一部分,20行。

x <- [7.6,2.2,1.1,4.7,8.6,7.5,7.5,29.9,5.0,3.0,2.4,1.5,14.9,3.9,3.7,3.2,5.0,1.7,2.9,2.3]

这里是功能说明

  1. 幂律:y=A*x^-(u)
  2. 指数:y=B*exp^(-βx)

现在,我想使用 MLE(最大似然法)得到 u 用于幂律,β 用于指数分布。

#set likelihood function of power law
pl <- function(u){-n*log(u-1)-n*(u-1)*log(min(x))+u*sum(log(x))}

#set likelihood function of exponential distribution
ex <- function(β){-n*log(β)+β*sum(x)}

这些功能对吗?

使用mle2()获取参数:

#get the parameter u of power law
s1 <- mle2(pl,start = list(u=2),data = list(x))
summary(s1)
#get the parameter lamda of exponential distribution
s2 <- mle2(ex,start = list(β=2),data = list(x))
summary(s2)

现在有两个问题:

  1. 如何确定哪一个最适合模型

    使用 confint() 可以获得 95% CI,如何获得两个模型的 Rsquared 和 AIC(Akaike weigths)?

  2. 在得到最适合数据的拟合后,如何在原始数据上方绘制拟合图?

我在 Windows 7 中使用 R.3.2.2。

【问题讨论】:

    标签: r exponential mle power-law


    【解决方案1】:

    几乎如您所料。您尚未指定数据的条件分布,因此我将假设正态性。 (鉴于此,您还可以使用nls() -- 最小二乘 最大似然估计来进行正常的同方差响应),尽管mle2 为使用优化器等提供了更多的空间。 )

    我打算用公式界面,如果你的模型不太复杂的话,这个界面很方便。

    幂律

    mle2(y~dnorm(mean=A*x^(-mu),sd=exp(logsd),
         start=list(A=...,mu=...,logsd=...),
         ## no need for list() if mydata is already data.frame
         data=mydata)   
    

    指数

      mle2(y~dnorm(mean=B*exp(-beta*x),sd=exp(logsd),
         start=list(B=...,beta=...,logsd=...),
         data=mydata)  
    

    ...start 中的元素是任何 合理 的起始值。鉴于您上面的数据,这些方法应该在数据的子集上合理地工作。但是,它们可能在 1000 万次观测中表现不佳。我会考虑使用

    glm(y~x,family=gaussian(link="log"),data=mydata)
    

    拟合指数曲线和

    glm(y~log(x),family=gaussian(link="log"),data=mydata)
    

    拟合幂律曲线。

    【讨论】:

    • x=[0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9]
    • x=[0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9],我的数据是这样的,但mydata的长度是1000万.还有一个问题,A、mu 和 logsd 可以是我在开始列表中输入的任何数字吗?
    • 你可能认为 mydata 是 data.frame 有 2 列名为 x 和 y?
    • 能否编辑您的问题以包含此信息?
    • 对不起,我对mydata的描述不好,我在问题描述中上传了一部分mydata,请再看一遍,谢谢!
    猜你喜欢
    • 1970-01-01
    • 2021-06-23
    • 1970-01-01
    • 1970-01-01
    • 2015-01-24
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多