﻿rm(list = ls())  # ワークスペースをクリア
hm <- function(m, x) { (x/m)^(m/2) * exp((m - x)/2) }  # 関数hm(x)
# 0<lambda<1に対してhm(x)=lambdaをみたすx(<m)を求める関数
hm_inv_left <- function(m, lambda) {
  xmin <- 0 ; xmax <- m; 
  xmid <- (xmin + xmax) / 2
  for (i in 1:50) {  # 二分法：開始
    if (hm(m, xmid) < lambda) {
      xmin <- xmid
    } else {
      xmax <- xmid
    }
    xmid <- (xmin + xmax) / 2
  }  # 二分法：終了
  return(xmid)
}
# 0<lambda<1に対してhm(x)=lambdaをみたすx(>m)を求める関数
hm_inv_right <- function(m, lambda) {
  xmin <- m ; xmax <- 2 * m
  # hm(x)=lambdaをみたすxが含まれる区間[k*m,(k+1)*m]を探索
  while (hm(m, xmax) > lambda) { 
    xmin <- xmax ; xmax <- xmax + m
  }
  xmid <- (xmin + xmax) / 2
  for (i in 1:50) {  # 二分法：開始
    if (hm(m, xmid) < lambda) {
      xmax <- xmid
    } else {
      xmin <- xmid
    }
    xmid <- (xmin + xmax) / 2
  }  # 二分法：終了   
  return(xmid)
}
#----------------------------------------------------------------
# αとmに対し，hm(xL)=hm(xR)，かつ
# ∫[0,xL] f^(m)_s0^2(x) dx + ∫[xR,∞] f^(m)_s0^2(x) dx = α
# をみたすxL = hmL_inv_lambda, xR = hmR_inv_lambdaを計算し
# 一覧表を出力する (m = 1, 2, ..., m_max)
#----------------------------------------------------------------
alpha <- 0.05  # 有意水準
m_max <- 30  # m = 1, 2, ..., m_max.
m_vec <- numeric(m_max)
xleft <- numeric(m_max)
xright <- numeric(m_max)
for(m in 1:m_max) {
  lambdamin <- 0 ; lambdamax <- 1
  for(i in 1:50) {  # 二分法：開始
    lambdamid <- (lambdamin + lambdamax) / 2
    xL <- hm_inv_left(m, lambdamid)
    xR <- hm_inv_right(m, lambdamid)
    # IL = ∫[0,xL] f^(m)_s0^2(x) dx, IR = ∫[xR,∞] f^(m)_s0^2(x) dx
    IL <- pchisq(xL, df = m)
    IR <- pchisq(xR, df = m, lower.tail = FALSE)
    if(IL + IR > alpha) {
      lambdamax <- lambdamid
    } else {
      lambdamin <- lambdamid
    }
  }  # 二分法：終了
  m_vec[m] <- m
  xleft[m] <- xL ; xright[m] <- xR
}
data.frame(m = m_vec, hmL_inv_lambda = xleft, hmR_inv_lambda = xright)
#------------------------------------------------------------------
# 参考：hm(xL)=hm(xR)=λをみたすxL<m<xRに対し
# I(λ) = ∫[0,xL] f^(m)_s0^2(x) dx + ∫[xR,∞] f^(m)_s0^2(x) dx
# のグラフを描画
#------------------------------------------------------------------
m <- 3  # m = n - 1
N <- 100  # [0,1]の分割数 (0<=λ<=1)
IL_lambda <- 0 ; IR_lambda <- 0 ; I_lambda <- 0
for(i in 1:N) {  # forループ
  lambda <- i / N
  xL <- hm_inv_left(m, lambda)
  xR <- hm_inv_right(m, lambda)
  # I_L = ∫[0,xL] f^(m)_s0^2(x) dx, I_R = ∫[xR,∞] f^(m)_s0^2(x) dx
  I_L <- pchisq(xL, df = m)  
  I_R <- pchisq(xR, df = m, lower.tail = FALSE)
  IL_lambda <- c(IL_lambda, I_L) ; IR_lambda <- c(IR_lambda, I_R)
  I_lambda <- c(I_lambda, I_L + I_R)
}
# 積分の値を順にプロットし，直線で結ぶ
x <- seq(0, 1, by = 1/N)  # x=(0, 1/N, 2/N, ..., 1)
plot(x = x, I_lambda, main = "m = 3", xlim = c(0, 1), ylim = c(0, 1),
     xlab = "lambda", ylab = "I(lambda)，IL(lambda)，IR(lambda)", 
     type = "l", lwd = 2.5, lty = 1, col = 4)
par(new = T)
plot(x = x, IL_lambda, main = "m = 3", xlim = c(0, 1), ylim = c(0, 1), 
     xlab = "lambda", ylab = "I(lambda)，IL(lambda)，IR(lambda)", 
     type = "l", lwd = 2.5, lty = 3, col = 1, ann = F)
par(new = T)
plot(x = x, IR_lambda, main = "m = 3", xlim = c(0, 1), ylim = c(0, 1), 
     xlab = "lambda", ylab = "I(lambda)，IL(lambda)，IR(lambda)",
     type = "l", lwd = 2.5, lty = 2, col = 1, ann = F)