# 対数正規確率紙に累積相対度数をプロットする
log.normal.probability.plot <- function(data=data.df, 
                                        log.buff=c(0.2, 0.2), 
                                        points.pch=NULL, 
                                        points.color=NULL,
                                        lines.color=NULL, 
                                        trimming=NULL)   # データベクトル
{
  x <- data[,1]
  x.condition <- data[,2]
  log.x.min <- min(log10(x))
  log.x.max <- max(log10(x))
  data.condition <- unique(x.condition)

  log.axis <- function(z)                   # 対数目盛りの軸を描く
  {
    z <- floor(log10(z))                    # 対数にしたときの整数部
    log.min <- min(z)                       # 最小値
    z2 <- 1:10*10^log.min                   # 一番左側に描かれる目盛り数値のセット
    n <- max(z)-log.min                     # 10 倍しながら順次，右の位置に目盛りを描く
    z2 <- rep(z2, n+1)*10^rep(0:n, each=10) # 対数目盛り位置の数値
    log.z2 <- log10(z2)                     # 目盛りを描く位置
    axis(1, at=log.z2, labels=z2)           # log.z2 の位置に，z2 という数値を描く
    abline(v=log.z2, col="gray")            # 垂直格子線を入れる
  }

  log.x <- log10(sort(x))
  log.x <- append(min(log.x)-log.buff[1], log.x)
  log.x <- append(log.x, max(log.x)+log.buff[2])
  probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, # 縦軸の目盛り
             60, 70, 80, 90, 95, 99, 99.9, 99.99)/100
  plot(log.x[c(1, length(log.x))], qnorm(probs[c(1,17)]),      # 枠組み
       type="n", xaxt="n", yaxt="n",
       xlab="Observed Value", ylab="Cumulative Percent",
       main="Log Normal Probability Paper")
  log.axis(x)                               # 横軸（対数目盛）を描く
  axis(2, qnorm(probs), probs*100)          # 縦軸（正規確率目盛）を描く
  abline(h=qnorm(probs), col="grey")        # 水平格子線を描く

  for( i in data.condition) {
    log.x <- log10(sort(x[x.condition==i]))
    n <- length(log.x)
    y <- 1:n/(n+1)
    k <- match(i, data.condition)
    if(!is.null(points.pch)&&!is.null(lines.color)) {
      points(log.x, qnorm(y), pch=points.pch[k], col=points.color[k])
    } else {
      points(log.x, qnorm(y), pch=20)
    }
    if( (!is.null(trimming)) & (length(trimming)==2) ) {
      i.min <- round(n*trimming[1])
      i.max <- round(n*trimming[2])
    } else {
      i.min <- 1
      i.max <- n
    }
    y.to.fit <- y[i.min:i.max]
    log.x.to.fit <- log.x[i.min:i.max]
    lm.fit <- lm(qnorm(y.to.fit)~log.x.to.fit)
    print(summary(lm.fit)$coefficients)
    print(summary(lm.fit)$r.squared)
    if(!is.null(lines.color)) {
      abline(lm.fit, col=lines.color[k])
    } else {
      abline(lm.fit)
    }
  }
}