9

私はRでいくつかのデータ分析を行っており、データを3パラメーターのワイブル分布に適合させる方法を見つけようとしています。2パラメーターのワイブルでそれを行う方法を見つけましたが、3パラメーターでそれを行う方法を見つけるのに不足しています。

パッケージのfitdistr関数を使用してデータを適合させる方法は次のとおりです。MASS

y <- fitdistr(x[[6]], 'weibull')

x[[6]]は私のデータのサブセットであり、yはフィッティングの結果を保存している場所です。

4

2 に答える 2

8

まず、FAdistパッケージを確認することをお勧めします。rweibull3しかし、それはからに行くのはそれほど難しいことではありませんrweibull：

> rweibull3
function (n, shape, scale = 1, thres = 0) 
thres + rweibull(n, shape, scale)
<environment: namespace:FAdist>

同様にdweibull3からdweibull

> dweibull3
function (x, shape, scale = 1, thres = 0, log = FALSE) 
dweibull(x - thres, shape, scale, log)
<environment: namespace:FAdist>

だから私たちはこれを持っています

> x <- rweibull3(200, shape = 3, scale = 1, thres = 100)
> fitdistr(x, function(x, shape, scale, thres) 
       dweibull(x-thres, shape, scale), list(shape = 0.1, scale = 1, thres = 0))
      shape          scale          thres    
    2.42498383     0.85074556   100.12372297 
 (  0.26380861) (  0.07235804) (  0.06020083)

編集：コメントで述べたように、このようにディストリビューションを適合させようとすると、さまざまな警告が表示されます

Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573,  : 
  non-finite finite-difference value [3]
There were 20 warnings (use warnings() to see them)
Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573,  : 
  L-BFGS-B needs finite values of 'fn'
In dweibull(x, shape, scale, log) : NaNs produced

私にとっては、最初はそれだけNaNs producedでしたが、初めて見たわけではないので、見積もりが良かったのであまり意味がないと思いました。いくつか検索したところ、それは非常に人気のある問題のようで、原因も解決策も見つかりませんでした。stats4パッケージと機能を使用することもできますmle()が、問題もあるようです。しかし、私が数回チェックしたdanielmedicによるコードの修正バージョンを使用することを提案できます。

thres <- 60
x <- rweibull(200, 3, 1) + thres

EPS = sqrt(.Machine$double.eps) # "epsilon" for very small numbers

llik.weibull <- function(shape, scale, thres, x)
{ 
  sum(dweibull(x - thres, shape, scale, log=T))
}

thetahat.weibull <- function(x)
{ 
  if(any(x <= 0)) stop("x values must be positive")

  toptim <- function(theta) -llik.weibull(theta[1], theta[2], theta[3], x)

  mu = mean(log(x))
  sigma2 = var(log(x))
  shape.guess = 1.2 / sqrt(sigma2)
  scale.guess = exp(mu + (0.572 / shape.guess))
  thres.guess = 1

  res = nlminb(c(shape.guess, scale.guess, thres.guess), toptim, lower=EPS)

  c(shape=res$par[1], scale=res$par[2], thres=res$par[3])
}

thetahat.weibull(x)
    shape     scale     thres 
 3.325556  1.021171 59.975470
于 2012-08-05T16:20:31.893 に答える
1

別の方法：パッケージ「lmom」。Lモーメント法による推定

library(lmom)
thres <- 60
x <- rweibull(200, 3, 1) + thres
moments = samlmu(x, sort.data = TRUE)
log.moments <- samlmu( log(x), sort.data = TRUE )
weibull_3parml <- pelwei(moments)
weibull_3parml
zeta      beta     delta 
59.993075  1.015128  3.246453

しかし、このパッケージまたは上記のソリューションで適合度統計を行う方法がわかりません。適合度統計を簡単に実行できるその他のパッケージ。とにかく、ks.testやchisq.testのような代替手段を使用できます

于 2018-09-25T17:57:50.803 に答える