January 26, 2015

指数分布の確率プロット

2年ぐらい前にお勉強用に書いたコードを再利用するために発掘したので、共有しておきます。

“exp”の箇所を好きな分布に書きかえれば流用できるかと思います。

exponential.plot <- function(x){

  qq <- function(x){

    n <- length(x)

    qexp((1:n - 1/2)/n, rate=mean(x))[order(order(x))] # ミラーランク

  }

  

  probs<- c(0.001, 0.01, 0.1, 1, 5, 10, 20, 30, 50, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999) / 100

  

  plot(c(min(x), max(x)),  qexp(c(probs[1], probs[length(probs)])), type="n", xaxt="n", yaxt="n", 

       ylab="Cumulative Percent", xlab="Observed Value")

  abline(h=qexp(probs), v=as.integer(seq(min(x), max(x))), col="grey")

  abline(h=qexp(0.5))

  axis(side=2, at=qexp(probs[c(4,8:17)]), las=1, labels=as.character(probs[c(4,8:17)]*100), cex.axis=0.8)

  axis(side=1, as.integer(seq(min(x), max(x))), cex.axis=0.8)

  points(x, qq(x))

  

  model <- lm(qq(x) ~ x - 1)

  abline(model)

  print(paste("r-squared",summary(model)$r.squared,sep=":"))

  # Maximum likelihood estimation

  print(paste("rate",mean(x),sep=":"))

}



rexp.data <- rexp(100,rate=1)

exponential.plot(rexp.data)

出力結果は、

exponential_plot

となります。

参考

© gepuro 2013