【R言語】母分散の信頼区間の計算

R言語で母分散の信頼区間を算出する関数やその実行例について解説します。

実行例では、データセットから母分散の信頼区間を計算するだけでなく、母分散の信頼区間の意味を図で説明します。

この記事で扱うプログラミングコードは以下からダウンロードできます。

R言語　母分散の信頼区間

1 ファイル 1.20 KB

ダウンロード

母分散の信頼区間差については、次の記事を参照してください。

母分散の信頼区間

分散の信頼区間とR言語の関数についてまとめました。

母分散の信頼区間の概要について列挙しています。

母分散の信頼区間

正規母集団$N(\mu, \sigma^2)$からの大きさ$n$の無作為標本$x_1, \ldots, x_n$が得られたときの、母分散$\sigma^2$の$100(1- \alpha)$%信頼区間を紹介する。以降、標本平均を$\bar{x} = (1 /n) \sum_{i = 1}^n x_i$、不偏標本分散を$u^2 = \sum_{i=1}^n (x_i - \bar{x})^2 / (n-1)$、自由度$n$のカイ二乗分布の上側$\alpha$点を$\chi_{n, \alpha}^2$とする。

母分散の信頼区間

母分散$\sigma^2$の$100(1 - \alpha)$%信頼区間は次で与えられる。

\begin{align} \label{eq2} \left[ \cfrac{(n - 1)u^2}{\chi_{n-1, \alpha / 2}^2 }, \cfrac{(n-1) u^2}{\chi_{n-1, 1 -\alpha / 2}^2} \right]. \tag{2}\end{align}

関数varTest

パッケージEnvStatsの関数varTestを用いることで母分散の検定を行うことができます。

以下varTestとその引数です。

varTest(x, alternative = "two.sided", conf.level = 0.95, sigma.squared = 1, data.name = NULL)

関数varTestの引数
引数	説明
x	numeric型のベクトル（観測値）。NAやNaN、infにも適用できるが除外される。
alternative	character型。"two.sided"で両側検定、"greater"で右片側検定、"less"左片側検定を指定できる。
conf.level	信頼区間の信頼度数。
sigma.squared	帰無仮説のときの母分散の値。
data.name	character型。検定に用いるデータの名前。

実行例

母分散の信頼区間の実行例について紹介します。

母分散の信頼区間の算出のためにパッケージEnvStats、グラフの作成の際にパッケージggplot2、gridExtraを用います。あらかじめインストールしておいてください。

install.packages("ggplot2")
install.packages("EnvStats")
library(ggplot2)
library(EnvStats)

install.packages("ggplot2")

install.packages("EnvStats")

library(ggplot2)

library(EnvStats)

次の木に関するデータセットirisを例に母分散の信頼区間の計算を行っていきます。

> dataset <- iris
> head(dataset)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

> dataset <- iris

> head(dataset)

Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1 5.1 3.5 1.4 0.2 setosa

2 4.9 3.0 1.4 0.2 setosa

3 4.7 3.2 1.3 0.2 setosa

4 4.6 3.1 1.5 0.2 setosa

5 5.0 3.6 1.4 0.2 setosa

6 5.4 3.9 1.7 0.4 setosa

SpeciesがsetosaであるSepal.Lengthのヒストグラムは以下の通りです。

boxPlot <- ggplot(dataset[dataset$Species == "setosa", ], aes(x = 1, y = Sepal.Length)) +
  geom_boxplot() +
  scale_y_continuous() +
  coord_flip() +
  geom_point(position = position_dodge(width = 0.9), alpha = 0.3, size = 1)

histPlot <- ggplot(dataset[dataset$Species == "setosa", ], aes(x = Sepal.Length)) +
  geom_histogram(bins = 15) +
  scale_y_continuous() 

grid.arrange(histPlot, boxPlot)

boxPlot <- ggplot(dataset[dataset$Species == "setosa", ], aes(x = 1, y = Sepal.Length)) +

geom_boxplot() +

scale_y_continuous() +

coord_flip() +

geom_point(position = position_dodge(width = 0.9), alpha = 0.3, size = 1)

histPlot <- ggplot(dataset[dataset$Species == "setosa", ], aes(x = Sepal.Length)) +

geom_histogram(bins = 15) +

scale_y_continuous()

grid.arrange(histPlot, boxPlot)

SpeciesがsetosaであるSepal.Lengthのヒストグラム

SpeciesがsetosaであるSepal.Lengthの母分散の信頼区間の計算を行いたいと思います。

setosaに限定する理由は、すべてのSpeciesを含めてしまうと正規性が成り立たなくなるためです。ヒストグラムを描けば分かりますが、全てを含めた場合Sepal.Lengthのヒストグラムは多峰性を持ちます。

SpeciesがsetosaであるSepal.Lengthの95%信頼区間を算出するには、次のように関数varTestの引数conf.levelを0.95を与えて後ろに$conf.intを付けます。

conf_level <- 0.95

ci <- varTest(dataset$Sepal.Length[dataset$Species == "setosa"], conf.level = conf_level)$conf.int

conf_level <- 0.95

ci <- varTest(dataset$Sepal.Length[dataset$Species == "setosa"], conf.level = conf_level)$conf.int

ciを参照すると信頼区間の下限と上限、また信頼水準が格納されているのが分かります。

> ci
       LCL        UCL 
0.08669881 0.19293982 
attr(,"conf.level")
[1] 0.95

> ci

LCL UCL

0.08669881 0.19293982

attr(,"conf.level")

[1] 0.95

信頼区間の上限と下限を取得するには、次のようにciの要素を名前で取得すれば大丈夫です。

ci_lower <- ci["LCL"] #信頼区間の下限
ci_upper <- ci["UCL"] #信頼区間の上限

1 2	ci_lower <- ci["LCL"] #信頼区間の下限 ci_upper <- ci["UCL"] #信頼区間の上限

計算した信頼区間を他のファイルに保存したいときは、次のように関数data.frameでデータフレームに計算した値を格納しwrite.csvなどで他のファイルへ保存します。

resultTable <- data.frame(sample.variance = var(dataset$Sepal.Length[dataset$Species == "setosa"]), 
                          c.i.lower = ci["LCL"], c.i.upper = ci["UCL"], row.names = NULL)
write.csv(resultTable, "母分散の信頼区間.csv", row.names = FALSE)

resultTable <- data.frame(sample.variance = var(dataset$Sepal.Length[dataset$Species == "setosa"]),

c.i.lower = ci["LCL"], c.i.upper = ci["UCL"], row.names = NULL)

write.csv(resultTable, "母分散の信頼区間.csv", row.names = FALSE)

母分散の信頼区間

上で算出した信頼区間をグラフにすると下の画像のようになります。

dataset_ci <- data.frame(x = c(mean(dataset$Sepal.Length[dataset$Species == "setosa"]) +
                                 c(sd(dataset$Sepal.Length[dataset$Species == "setosa"]), 
                                   -sd(dataset$Sepal.Length[dataset$Species == "setosa"])),
                               mean(dataset$Sepal.Length[dataset$Species == "setosa"]) +
                                 c(sqrt(ci), -sqrt(ci))),
                         name = factor(c(rep("SD", 2), paste0(100 * conf_level, paste0("% c.i.", rep(c("upper", "lower"), 2)))),
                                       levels = c("SD", paste0(100 * conf_level, "% c.i.", c("upper", "lower")))))

ggplot(dataset[dataset$Species == "setosa", ], aes(x = 1, y = Sepal.Length)) + 
  geom_dotplot(binaxis = "y", stackdir = "center") +
  stat_summary(aes(x = 1, y = Sepal.Length), 
               fun.min = function(x)  mean(x) - sd(x), fun.max = function(x) mean(x) + sd(x), 
               geom = 'errorbar', width = 0.15, color = "#00BA38", size = 1) +
  stat_summary(aes(x = 0.95, y = Sepal.Length), 
               fun.min = function(x)  mean(x) - sqrt(ci["LCL"]), fun.max = function(x) mean(x) + sqrt(ci["LCL"]), 
               geom = 'errorbar', width = 0.1, color = "#619CFF", size = 1) + 
  stat_summary(aes(x = 1.05, y = Sepal.Length), 
               fun.min = function(x)  mean(x) - sqrt(ci["UCL"]), fun.max = function(x) mean(x) + sqrt(ci["UCL"]), 
               geom = 'errorbar', width = 0.1, color = "#F8766D", size = 1) +
  geom_text(data = dataset_ci,
            aes(x = 0.7, y = x, color = name, label = paste0(c(rep("SD", 2), rep(c("c.i.lower", "c.i.upper"), 2)), ": ", round(x, 3)))) +
  labs(color = "") +
  scale_colour_manual(values = c("#00BA38", "#619CFF", "#F8766D"))

dataset_ci <- data.frame(x = c(mean(dataset$Sepal.Length[dataset$Species == "setosa"]) +

c(sd(dataset$Sepal.Length[dataset$Species == "setosa"]),

-sd(dataset$Sepal.Length[dataset$Species == "setosa"])),

mean(dataset$Sepal.Length[dataset$Species == "setosa"]) +

c(sqrt(ci), -sqrt(ci))),

name = factor(c(rep("SD", 2), paste0(100 * conf_level, paste0("% c.i.", rep(c("upper", "lower"), 2)))),

levels = c("SD", paste0(100 * conf_level, "% c.i.", c("upper", "lower")))))

ggplot(dataset[dataset$Species == "setosa", ], aes(x = 1, y = Sepal.Length)) +

geom_dotplot(binaxis = "y", stackdir = "center") +

stat_summary(aes(x = 1, y = Sepal.Length),

fun.min = function(x) mean(x) - sd(x), fun.max = function(x) mean(x) + sd(x),

geom = 'errorbar', width = 0.15, color = "#00BA38", size = 1) +

stat_summary(aes(x = 0.95, y = Sepal.Length),

fun.min = function(x) mean(x) - sqrt(ci["LCL"]), fun.max = function(x) mean(x) + sqrt(ci["LCL"]),

geom = 'errorbar', width = 0.1, color = "#619CFF", size = 1) +

stat_summary(aes(x = 1.05, y = Sepal.Length),

fun.min = function(x) mean(x) - sqrt(ci["UCL"]), fun.max = function(x) mean(x) + sqrt(ci["UCL"]),

geom = 'errorbar', width = 0.1, color = "#F8766D", size = 1) +

geom_text(data = dataset_ci,

aes(x = 0.7, y = x, color = name, label = paste0(c(rep("SD", 2), rep(c("c.i.lower", "c.i.upper"), 2)), ": ", round(x, 3)))) +

labs(color = "") +

scale_colour_manual(values = c("#00BA38", "#619CFF", "#F8766D"))

母分散の信頼区間の意味

緑線はsetosa種のSepal.Lengthの標本平均±標準偏差
赤線はsetosa種のSepal.Lengthの標本平均±標準偏差の信頼区間の上限
青線はsetosa種のSepal.Lengthの標本平均±標準偏差の信頼区間の下限

を表しています。

母分散の信頼区間がだいたい0.08から0.19であることから、その平方根である0.29から0.44が標準偏差の信頼区間となります。

データのばらつきの指標である標本平均±標準偏差は上の青から赤の範囲にあることが分かります。（標本平均は考慮していません。）

まとめ

R言語で母分散の信頼区間を計算する関数やその実行例をみていきました。

パッケージEnvStatsの関数varTestを用いることで母分散の信頼区間を計算することができます。

usagi-san

統計学とゲームとかをメインに解説していくよ。数式とかプログラミングコードにミスがあったり質問があったりする場合はコメントで受け付けます。すぐに対応します。

2024/04/28

【R言語】F検定テンプレートスクリプト

2024/04/28

【R言語】経験密度関数・経験分布関数のプロット　関数densityとecdfの使い方

2024/04/10

【R言語】ベータ関数とガンマ関数　関数beta, gammaの使い方

usagi-sanの記事をもっと見る

-R言語
-R言語, プログラミング, 統計学, 統計解析

comment コメントをキャンセル

: R言語
【R言語】関数t.testを用いた2標本t検定・母平均の差の検定
母平均の差についての検定を行う方法を解説していきます。 R言語を用いることで、母平均の差の検定である2標本t検定が簡単に実行できます。 1標本の時と同様に、検定結果などは自分で計算する必要がなく非常に ...

: R言語
【R言語】加重平均の計算　重み付きの平均　記述統計
R言語で重み付きの平均である加重平均を計算する関数やその使用例について紹介していきます。実行例では、加重平均の具体的な計算例について見ていき、どのような場面で用いられるのか理解を深めていきます。こ ...

: R言語
【R言語】リッジ回帰・ラッソ回帰　glmnetを用いたL1, L2正則化
R言語でリッジ回帰・ラッソ回帰を実行する関数やその実行例を紹介します。この記事では、モデルの過適合を回避するための正則化（regularization）をRで行い、実際にテストデータに対する予測精度 ...