2015-01-05

確率プロットの作成

R によるグラフ作成の復習の第三弾の今回は、確率プロットです。以前の記事は下記の通りです。

今回も、ファイル選択ダイアログを加えました。例に使用したデータはボックスプロットで使用したデータと同じで、sample_box.csv.zip からダウンロードできるようにしてあります。ご興味があれば、展開してご利用下さい。

require("tcltk")

###############################################################################
# chart.cdf
###############################################################################
chart.cdf <- function(tbl, x.val = "RAW_VALUE", col.by = "condition",
                      gTitle = "", x.label = "measured data") {
  x <- tbl[, x.val]
  x <- x[!is.na(x)]
  n <- length(x)
  x <- sort(x)
  y <- (1:n - 0.375) / (n + 0.25)
  
  probs <- c(0.001, 0.01, 0.1, 1, 5, 10, 20, 30, 40, 50,
             60, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999) / 100
  m <-  length(probs)
  
  x_min <- x[1]
  x_max <- x[n]
  x_pad <- (x_max - x_min) * 0.4
  
  y_min <- qnorm(probs[1])
  y_max <- qnorm(probs[m])
  y_stp <- (y_max - y_min) * 0.1
  
  plot(c(x_min, x_max + x_pad), c(y_min, y_max),
       type = "n", yaxt = "n",
       xlab = x.label, ylab = "Cumulative Percent",
       main = gTitle, cex.main = 1)
  axis(2, qnorm(probs), probs * 100)
  abline(h = qnorm(probs), lty = 3, col = "#C0C0C0")
  abline(h = 0, lty = 1, col = "#808080")
  
  list.by <- sort(unique(tbl[, col.by]))
  k <- length(list.by)
  for (i in 1:k) {
    val.by <- list.by[i]
    x <- tbl[tbl[, col.by] == val.by, x.val]
    x <- x[!is.na(x)]
    n <- length(x)
    x <- sort(x)
    y <- (1:n - 0.375) / (n + 0.25)
    
    total <- length(x)
    x.median <- median(x)
    legend_label <- paste(val.by, ':', total, "/", x.median)
    
    points(x, qnorm(y), type='o', col = i, pch = i, cex = 0.7)
    
    legend_x <- x_max + x_pad * 0.05
    legend_y <- y_max - y_stp * (i - 1)
    legend_lty = 1
    legend(legend_x, legend_y, legend_label, col = i,
           lty = legend_lty, pch = i, ncol = 1, cex = 0.7, box.lty = 0)
  }
}

###############################################################################
# open.png
###############################################################################
open.png <- function(img.name, pic.width = 640, pic.height = 400) {
  png(filename = img.name, width = pic.width, height = pic.height,
      units = "px", pointsize = 14, bg = "white", res = NA)
}

###############################################################################
# MAIN
###############################################################################

# READ CSV DATA
file.name <- tclvalue(tkgetOpenFile(filetypes = "{{CSV} {.csv}}", title = "Open CSV file"))
tbl.data <- read.csv(file.name, header = T, as.is = T)

# generate CDF PLOT
head.y <- names(tbl.data)[1]
head.x <- names(tbl.data)[2]
title.box <- paste("SAMPLE (n = ", length(tbl.data[, head.x]), ")", sep = "")
label.x <- "Measurement"
chart.cdf(tbl.data, head.x, head.y, title.box, label.x)

# generate PNG FILE
open.png("sample_cdf.png")
chart.cdf(tbl.data, head.x, head.y, title.box, label.x)
dev.off()

# ---
# END PROGRAM

Fedora での実行例です。

$ R

R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"
...
(省略)
...
 'q()' と入力すれば R を終了します。 

> source("cdf.R")
 要求されたパッケージ tcltk をロード中です 
> 

スクリプトを実行するとファイル選択ダイアログが表示されます。

画面に表示される確率プロットです。凡例には、データ数とメディアンを表示しています。

カレントディレクトリに PNG 形式で保存された確率プロットです。

 

0 件のコメント: