1

ガンマ分布曲線を他のべき乗則のような曲線のプロットに重ね合わせる必要があります。まず、ヒストグラムのドットポイントを両対数スケールでプロットします

  plot(log(pp$mids),log(pp$density))

次に、外部関数gamma()を呼び出してガンマ分布曲線を重ね合わせたいと思います。

  gamma <- function(X)
  {
  n <- length(X)
  theta<-var(hh2$V1)/mean(hh2$V1)
  kappa<-mean(hh2$V1)/theta
  y<-rgamma(n,kappa,theta)
  xx<-hist(y,plot=F)
  curve(log(xx$density),add=T,col='violet',type='l')
  return( c(kappa) ) 
  } 

しかし、curve()はプロットするために真の曲線を必要とするため、これはエラーを返します。これどうやってするの?

4

1 に答える 1

1

コードの多少機能するバリアントは次のとおりです。

(おそらく)データが次のように構造化された例を生成します。

library(rmutil)  ## for rpareto
set.seed(101)
hh2 <- data.frame(V1=rpareto(1000, m=1, s=1.5))

初期ヒストグラム計算:

pp <- hist(hh2$V1,plot=FALSE)

gamma関数(組み込み関数をマスクするため、呼び出さない方がよい):

ghistfun <- function(x) {
    n <- length(x)
    scalepar <- var(x)/mean(x)
    shapepar <- mean(x)^2/var(x)
    y <- rgamma(n,shape=shapepar,scale=scalepar)
    xx <- hist(y,plot=FALSE)
    lines(log(xx$mids),log(xx$density),col="red")
    curve(dgamma(exp(x),shape=shapepar,scale=scalepar,log=TRUE),
        add=TRUE,col="blue")

    shapepar
}

n正確に同じサイズのデータ​​セットのランダムな変動を確認することに特に関心がない限り、データの長さだけを使用するのではなく、非常に大きな数を使用する方がよい場合があります。または、示されているように、を使用することもできます(最初は、の密度からの密度へcurve(dgamma(x,...))のスケーリングを許可する必要があると思いましたが、元の(ログに記録されていない)スケールでヒストグラムを計算してから、ビンの中点、あなたはする必要はありません...)xlog(x)

plot(log(pp$mids),log(pp$density))
ghistfun(hh2$V1)

ここに画像の説明を入力してください

于 2012-07-27T15:06:11.087 に答える