baseパケージの用意されている関数をR内で（条件付きで動作する範囲で）再現する試み。
自分の勉強用。
Rのt分布に関する組み込み関数を使ってt分布の確率密度関数を描出するとすると、
curve(dt(x,10),xlim = c(-4,4))
となる。
dt()の本体はcで書かれたソースが見つかるが、そういうのを完全再現するのではなくt分布の確率密度関数f(x) = Γ( (n+1)/2) / (√(n π) Γ(n/2) ) (1 + x2/n)^-( (n+1)/2) )から再現すると、
t.dist <- function(t, n){ gamma((n + 1) / 2) / (sqrt(n * pi) * gamma(n / 2)) * (1 + t^2 / n)^(-(n + 1) / 2) } curve(t.dist(x,10),xlim = c(-4,4))
で同じfig.1が描ける。
関数内で確率密度関数を再現しているだけ。
t検定でt値から無限大（t値がマイナスなら-無限大）までの積分をして両側検定でもとめるなら二倍することになっている。
普通はpt()を使って計算する。
例えば、単一サンプルのt検定（両側）を関数自作すると
# 単一サンプルのt検定 # muは比較対象となる数値 # xがデータ one.sample.t.test <- function(x, mu = 0){ # データ数 n <- length(x) # 自由度算出 df <- n - 1 # t値算出 t.value <- (mean(x) - mu) / (sd(x) / sqrt(n)) # エラー計算 err <- qt(0.975, df = n - 1) * sd(x) / sqrt(n) # 95%信頼区間の計算 mean.conf.upper <- mean(x) + err mean.conf.lower <- mean(x) - err # p値の計算 p.value <- 2 * pt(- abs(t.value), df) # 計算した値をリストにする ans <- list("t value" = t.value, 'df' = df, "p-value" = p.value, 'mean' = mean(x), "95 percent confidence interval" = c(mean.conf.lower, mean.conf.upper) ) # 結果を返す ans }
となり、実行すると
x <- c(0.06, 0.32, 0.28, -0.5, -0.22, -1.47, 0.79, 1.33,-0.02, -0.32) one.sample.t.test(x) $`t value` [1] 0.1046803 $df [1] 9 $`p-value` [1] 0.9189253 $mean [1] 0.025 $`95 percent confidence interval` [1] -0.5152537 0.5652537 # t.testで検算 t.test(x) One Sample t-test data: x t = 0.10468, df = 9, p-value = 0.9189 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: -0.5152537 0.5652537 sample estimates: mean of x 0.025
となる。
ちなみに、t.testみたいな形で自作関数の結果を返したかったらlistの名前と投入項目整理して、listのclassをhtestにするといけると思う。
んで、本題のpt()を使わずにp値を求めるには最初のt分布の確率密度関数を積分すればいいだけなので、
# 単一サンプルのt検定（両側） # muは比較対象となる数値 # xがデータ one.sample.t.test2 <- function(x, mu = 0){ # データ数 n <- length(x) # 自由度算出 df <- n - 1 # t値算出 t.value <- (mean(x) - mu) / (sd(x) / sqrt(n)) # エラー計算 err <- qt(0.975, df = n - 1) * sd(x) / sqrt(n) # 95%信頼区間の計算 mean.ci <- mean(x) + c(-err, err) # mean.conf.lower <- mean(x) - err # p値の計算 # fでt分布の確率密度関数を定義 f <- function(t, n = df){ gamma((n + 1) / 2) / (sqrt(n * pi) * gamma(n / 2)) * (1 + t^2 / n)^(-(n + 1) / 2) } # pがt値の絶対値から上側の積分値（片側確率を算出） p <- integrate(f, abs(t.value), Inf) # p*2でp値算出（積分の計算上微妙な誤差がある） p.value <- 2 * p$value # 計算した値をリストにする ans <- list("t value" = t.value, 'df' = df, "p-value" = p.value, 'mean' = mean(x), "95 percent confidence interval" = mean.ci ) # 結果を返す ans } one.sample.t.test2(x) $`t value` [1] 0.1046803 $df [1] 9 $`p-value` [1] 0.9189253 $mean [1] 0.025 $`95 percent confidence interval` [1] -0.5152537 0.5652537
でpt()使って計算した値と基本的にほぼ同じ値を返す。ほぼ同じ値というのは、integrate()の計算で誤差が微妙にでるから。
積分するときにt値の絶対値をとらないと、t値が0以上か未満でif計算しなければいけなくなるが、絶対値をとって+無限大で積分を計算すると簡単になる。確率密度関数の形がわかっていれば他の検定でもできるはず。
こうやってみて本当に積分した面積でp値出るんだなぁと実感した。手計算するときにはt分布表をつかって計算するので、疑うわけじゃないけど本当にこの面積がp値なのかややブラックボックス感が自分の中ではあった。
t分布のRの既存の関数を使わないでt検定をやろうと結構いろいろ考えたけど、信頼区間を計算するために使うqt()のうまいやり方が思いつかずqt()だけ残ってしまった。
cでかかれたコードをみてみたけど何でかよくわからない処理で、うまくR内で再現できない。
