fitdist 截断法线答案

【问题标题】：fitdist for truncated normalfitdist 截断法线
【发布时间】：2012-08-23 18:11:03
【问题描述】：

截断的法线由下式给出：

dtnorm<- function(x, mean, sd, a, b) {
dnorm(x, mean, sd)/(pnorm(b, mean, sd)-pnorm(a, mean, sd))
}
ptnorm <- function(x, mean, sd, a, b) {
(pnorm(x,mean,sd) - pnorm(a,mean,sd)) / 
  (pnorm(b,mean,sd) - pnorm(a,mean,sd))
}

拟合由下式给出：

fitdist( data, tnorm, method="mle",
                    start=list(mean=mapply("[[", results[1], 1),
                               sd=mapply("[[", results[1], 2)),
                    fix.arg=list(a=minLoose,b=maxLoose))

其中 results[i] 是一个矩阵，其中 fitdist 的 mle 结果使用 normal 而不是 tnormal。

我得到以下关于 tnorm 的结果：

mean=-0.00844725266454969, sd=0.012540928272073

而规范：

mean=0.00748402597402597, sd=0.00614293813955003

数据都大于 0 且小于 0.04，因此为 tnorm 获得的 mle 似乎不正确....有什么建议吗？

谢谢！

【问题讨论】：

标签： r

【解决方案1】：

您的数据都高于正常值（呃，而不是高于 0）这一事实与截断分布的最佳拟合“平均值”是否超过 0 几乎没有关系。您正在拟合数据的正态分布。截断的估计位置参数并不是真正的平均值，而是平均值在未经审查的数据集中的位置，其右尾的密度“形状”与您的数据相同。（这实际上是一个统计问题，而不是 R 问题。）

您可以在维基百科文章的时刻部分找到计算双重截断法线的期望值的公式： http://en.wikipedia.org/wiki/Truncated_normal_distribution 它很容易转化为对pnorm 和qnorm 的调用。

进一步思考：查看在包中处理截断分布的工具：“gamlss”和“gamlss.tr”。

【讨论】：

我同意。但是如果我想知道截断范围内 tnorm 分布的平均值，那么我可以将它与范数的平均值进行比较......对吗？
截断法线的位置参数与其均值不同。如果您想要样本均值，您只需使用均值（数据）。如果您想要具有特定估计参数的分布的理论平均值，那么您应该计算 x*ptnorm() 跨越可接受域的限制。既然你还真的没有完全定义问题，也就不多说了。
好的，谢谢，我知道了。但是我只能使用矩的方法计算可接受域内的 tnorm 的矩吗？我也想使用 mle 计算它们以获得不确定性，不仅是 tnorm 的时刻，还有可接受域内的 tnorm。
我不确定我是否理解您的最后评论，但也许我的答案的补充是对它的回应？您确实有 mu、sigma 的估计值，并且可能已经预先指定了 a 和 b。
所以我首先运行 fitdist for norm 以获取我的数据的法线的 mle 估计值，然后将这些参数提供给 mle 拟合到截断的法线。结果，我得到了 tnorm 的均值和 sd 的 mle 估计值，它最适合我在可接受域上的数据。下一个问题是除了矩量法（您建议）之外，我还想再次使用 mle 来获得 tnorm 的均值和 sd 的估计值，现在仅在可接受的域上，因为 mle 也给了我不确定性估计值。

【解决方案2】：

您可以使用此脚本的一部分来估计参数

rm(list=ls(all=TRUE))


dtnorm<- function(x, mean, sd, a, b) {
dnorm(x, mean, sd)/(pnorm(b, mean, sd)-pnorm(a, mean, sd))
}


simuls=5
simul_mat=matrix(nrow=simuls,ncol=6)
for(simul in 1:simuls) {
acm=rnorm(1)
acsd=runif(1)*2+0.5
limits=sort(acm+rnorm(2))

all=limits[1]
aul=limits[2]



x=rnorm(10000)*acsd+acm
x=subset(x,x>all & x<aul)




    norm_parms<-function(parms){
    mp=parms[1]
    sdp=parms[2]^2
    ll=median(x)-parms[3]^2
    ul=median(x)+parms[4]^2

    xs=subset(x,x>ll & x<ul)
    ds=dtnorm(xs,mp,sdp,ll,ul)

    if(length(x)>5){
    do=rep(dnorm(-6),length(x)-length(xs))
    ds=c(ds,do)
    }
    if(length(x)<=5){
    ds=rep(dnorm(-9),length(x))
    }


    mll=-sum(log(ds))
    return(mll)
    }



bestv=Inf
methodss=c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN")
for(method in methodss){
try(bestc<-optim(par=c(0,1,1,1),norm_parms,method=method))
if(bestc$value<bestv) {best=bestc;bestv=bestc$value}
}



parms=best$par
mp=parms[1]
sdp=parms[2]^2
ll=median(x)-parms[3]^2
ul=median(x)+parms[4]^2
print(c(acm,acsd,all,aul))
print(c(mp,sdp,ll,ul))
print(best$value)
acparms=c(acm,acsd,sqrt(median(x)-all),sqrt(aul-median(x)))
acv=norm_parms(acparms) 
cnames=c("Actual a","Estimated a","Actual b","Estimated b","Actual optim","Best optim`")
simul_mat[simul,]=c(all,ll,aul,ul,best$value,acv)

cnames=c("Actual a","Estimated a","Actual b","Estimated b","Actual optim","Best optim`")
colnames(simul_mat)=cnames
print(simul_mat)

}

【讨论】：