2

切り捨てられた法線は次の式で与えられます。

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] は、tnormal の代わりに normal を使用した fitdist の mle 結果を含む行列です。

tnorm に対して次の結果が得られます。

mean=-0.00844725266454969, sd=0.012540928272073

一方、ノルムでは：

mean=0.00748402597402597, sd=0.00614293813955003

データはすべて 0 より大きく、0.04 より小さいため、tnorm で得られた mle は正しくないようです.... 何かアドバイスはありますか?

ありがとう！

4

2 に答える 2

5

データがすべて正規を上回っている (むしろ 0 を上回っている) という事実は、切り捨てられた分布への最適な適合の「平均」が 0 を超えるかどうかにはほとんど影響しません。正規分布の右裾を当てはめています。あなたのデータに。切り捨てられた場所の推定位置パラメーターは、実際には平均値ではなく、データと同じ密度の「形状」の右裾を持つ打ち切られていないデータセット内の平均値です。(これは実際には R の質問ではなく、統計に関する質問です。)

ウィキペディアの記事の瞬間セクションで二重に切り捨てられた法線の期待値を計算する式を見つけることができます : http://en.wikipedia.org/wiki/Truncated_normal_distributionpnormqnorm

さらに考えてみましょう: パッケージ内の切り捨てられたディストリビューションを操作するための機能を確認してください: 'gamlss' および 'gamlss.tr'。

于 2012-08-23T18:22:54.173 に答える
1

パラメータを推定するこのスクリプトの一部を使用できます

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)

}
于 2015-09-24T21:36:56.837 に答える