Python、Rで基本的な統計的推測まとめ

備忘録として、Python、Rでの基本的な統計的推測法のまとめを記します。

2標本のF検定

・2標本データから、それぞれの母分散に統計的に有意な差があるかどうかを調べたい時に利用する検定

Python

GitHub: https://github.com/Gin04gh/samples_py/blob/master/TwoSampleF-test.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
d1 = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
d2 = data[data["target"] == 2.0]["petal length (cm)"] # virginicaのPetal.Length
display(d1.head())
display(d2.head())
#   sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)    target
# 0 5.1 3.5 1.4 0.2 0.0
# 1 4.9 3.0 1.4 0.2 0.0
# 2 4.7 3.2 1.3 0.2 0.0
# 3 4.6 3.1 1.5 0.2 0.0
# 4 5.0 3.6 1.4 0.2 0.0
# 0    1.4
# 1    1.4
# 2    1.3
# 3    1.5
# 4    1.4
# Name: petal length (cm), dtype: float64
# 100    6.0
# 101    5.1
# 102    5.9
# 103    5.6
# 104    5.8
# Name: petal length (cm), dtype: float64

f = np.var(d1.values) / np.var(d2.values)
p = stats.f.cdf(f, len(d1)-1,  len(d2)-1)
print(p)
# 9.04487828439e-14

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/TwoSampleF-test.ipynb

# irisデータセット
head(iris)
d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
d2 <- subset(iris, Species=="virginica")$Petal.Length # virginicaのPetal.Length
head(d1)
head(d2)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7
# 6
# 5.1
# 5.9
# 5.6
# 5.8
# 6.6

var.test(d1, d2) # 等分散性の検定
#
#   F test to compare two variances
#
# data:  d1 and d2
# F = 0.099016, num df = 49, denom df = 49, p-value = 1.875e-13
# alternative hypothesis: true ratio of variances is not equal to 1
# 95 percent confidence interval:
#  0.05618945 0.17448557
# sample estimates:
# ratio of variances
#          0.0990164

ウィルコクソン（Wilcoxon）の順位和検定

・ノンパラメトリックな手法
・2標本データから、それぞれの母平均に統計的に有意な差があるかどうかを調べたい時に利用する検定
・ノンパラメトリック検定であるため、正規性を仮定できない場合に有効
・順位にのみ注目するため、極端に小さい値は極端に大きい値（外れ値）を含んでいる場合にも強い手法

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/TwoSampleWilcoxonTest.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
d1 = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
d2 = data[data["target"] == 2.0]["petal length (cm)"] # virginicaのPetal.Length
display(d1.head())
display(d2.head())
"""
    sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)    target
0   5.1 3.5 1.4 0.2 0.0
1   4.9 3.0 1.4 0.2 0.0
2   4.7 3.2 1.3 0.2 0.0
3   4.6 3.1 1.5 0.2 0.0
4   5.0 3.6 1.4 0.2 0.0
0    1.4
1    1.4
2    1.3
3    1.5
4    1.4
Name: petal length (cm), dtype: float64
100    6.0
101    5.1
102    5.9
103    5.6
104    5.8
Name: petal length (cm), dtype: float64
"""

p = stats.wilcoxon(d1, d2).pvalue # Wilcoxonの順位和検定（Mann-Whitney検定）
print(p) # p-value
"""
7.47869297957e-10
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/TwoSampleWilcoxonTest.ipynb

# irisデータセット
head(iris)
d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
d2 <- subset(iris, Species=="virginica")$Petal.Length # virginicaのPetal.Length
head(d1)
head(d2)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7
# 6
# 5.1
# 5.9
# 5.6
# 5.8
# 6.6

wilcox.test(d1, d2) # Wilcoxonの順位和検定（Mann-Whitney検定）
#
#   Wilcoxon rank sum test with continuity correction
#
# data:  d1 and d2
# W = 0, p-value < 2.2e-16
# alternative hypothesis: true location shift is not equal to 0

ウェルチ（Welch）のt検定

・2標本データから、それぞれの母平均に統計的に有意な差があるかどうかを調べたい時に利用する検定
・母分散が等しい（等分散性）と仮定できる場合は、スチューデントのt検定を利用するが、それが仮定できない場合にこちらのウェルチのt検定を利用する

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/TwoSampleWelchT-test.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
d1 = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
d2 = data[data["target"] == 2.0]["petal length (cm)"] # virginicaのPetal.Length
display(d1.head())
display(d2.head())
"""
sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)    target
0   5.1 3.5 1.4 0.2 0.0
1   4.9 3.0 1.4 0.2 0.0
2   4.7 3.2 1.3 0.2 0.0
3   4.6 3.1 1.5 0.2 0.0
4   5.0 3.6 1.4 0.2 0.0
0    1.4
1    1.4
2    1.3
3    1.5
4    1.4
Name: petal length (cm), dtype: float64
100    6.0
101    5.1
102    5.9
103    5.6
104    5.8
Name: petal length (cm), dtype: float64
"""

p = stats.ttest_ind(d1, d2, equal_var=False).pvalue # Welchのt検定（等分散性不明）
print(p) # p-value
"""
9.7138670617e-50
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/TwoSampleWelchT-test.ipynb

# irisデータセット
head(iris)
d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
d2 <- subset(iris, Species=="virginica")$Petal.Length # virginicaのPetal.Length
head(d1)
head(d2)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7
# 6
# 5.1
# 5.9
# 5.6
# 5.8
# 6.6

t.test(d1, d2, var.equal=FALSE) # Welchのt検定（等分散性不明）
#
#   Welch Two Sample t-test
#
# data:  d1 and d2
# t = -49.986, df = 58.609, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#  -4.253749 -3.926251
# sample estimates:
# mean of x mean of y
#     1.462     5.552

スチューデント（Student）のt検定

・2標本データから、母平均に統計的に有意な差があるかどうかを調べたい時に利用する検定
・母分散に等分散性が仮定できない場合には、ウェルチのt検定を利用する

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/TwoSampleStudentT-test.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
d1 = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
d2 = data[data["target"] == 2.0]["petal length (cm)"] # virginicaのPetal.Length
display(d1.head())
display(d2.head())
"""
sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)    target
0   5.1 3.5 1.4 0.2 0.0
1   4.9 3.0 1.4 0.2 0.0
2   4.7 3.2 1.3 0.2 0.0
3   4.6 3.1 1.5 0.2 0.0
4   5.0 3.6 1.4 0.2 0.0
0    1.4
1    1.4
2    1.3
3    1.5
4    1.4
Name: petal length (cm), dtype: float64
100    6.0
101    5.1
102    5.9
103    5.6
104    5.8
Name: petal length (cm), dtype: float64
"""

p = stats.ttest_ind(d1, d2).pvalue # Stundentのt検定（等分散性あり）
print(p) # p-value
"""
1.56412241589e-71
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/TwoSampleStudentT-test.ipynb

# irisデータセット
head(iris)
d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
d2 <- subset(iris, Species=="virginica")$Petal.Length # virginicaのPetal.Length
head(d1)
head(d2)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7
# 6
# 5.1
# 5.9
# 5.6
# 5.8
# 6.6

t.test(d1, d2, var.equal=TRUE) # Stundentのt検定（等分散性あり）
#
#  Two Sample t-test
#
# data:  d1 and d2
# t = -49.986, df = 98, p-value < 2.2e-16
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#  -4.252374 -3.927626
# sample estimates:
# mean of x mean of y
#     1.462     5.552

カイ二乗検定

・独立性の検定、一様性の検定、適合度検定など
・独立性の検定：観測されたデータを水準ごとに分けた時に、その水準間において関係があるかどうかを調べる
・適合度検定：観測されたデータを水準ごとに分けた時に、その水準ごとに仮定される確率（分布、割合）通りに観測値がばらついているかどうかを調べる
・一様性の検定：観測されたデータを水準ごとに分けた時に、その水準間において観測値のばらつきに偏りがない（一様に分布されている）かどうかを調べる
・期待度数が1以下のものが1つでもあったり、5以下のものが20%あるような場合はフィッシャーの正確確率検定を用いる

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/Chi-squaredTest.ipynb

import pandas as pd
import numpy as np
from scipy import stats

# データ
data = np.array([
    [43, 20],
    [13, 19]
])
"""
        治癒した 治癒しなかった
新薬 43　  20
偽薬 13   19
"""

squared, p, dof, ef = stats.chi2_contingency(data) # 独立性の検定（カイ二乗検定）
print("統計量: ", squared)
print("p-value: ", p)
print("自由度", dof)
print("期待度数")
print(ef)
"""
統計量:  5.60104188243
p-value:  0.017949799846
自由度 1
期待度数
[[ 37.13684211  25.86315789]
 [ 18.86315789  13.13684211]]
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/Chi-squaredTest.ipynb

# データ
data <- matrix(c(43, 20, 13, 19), ncol=2, nrow=2)
#        治癒した 治癒しなかった
# 新薬 43　  20
# 偽薬 13   19

chisq.test(data) # 独立性の検定（カイ二乗検定）
#
#  Pearson's Chi-squared test with Yates' continuity correction
#
# data:  data
# X-squared = 5.601, df = 1, p-value = 0.01795

コルモゴロフ・スミルノフ（Kolmogorov-Smirnov）検定

・サンプルデータから、母集団に正規性があるかどうかを調べたい時に利用する検定
・同じく母集団の正規性を調べる手法として、シャピロ・ウィルク検定がある

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/Kolmogorov-SmirnovTest.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
data = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
display(data.head())
"""
    sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)    target
0   5.1 3.5 1.4 0.2 0.0
1   4.9 3.0 1.4 0.2 0.0
2   4.7 3.2 1.3 0.2 0.0
3   4.6 3.1 1.5 0.2 0.0
4   5.0 3.6 1.4 0.2 0.0
0    1.4
1    1.4
2    1.3
3    1.5
4    1.4
Name: petal length (cm), dtype: float64
"""

data_z = stats.zscore(data.values) # 標準化
p = stats.kstest(data_z, "norm").pvalue # Kolmogorov-Smirnovの正規性の検定
print(p) # p-value
"""
0.15267909482
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/Kolmogorov-SmirnovTest.ipynb

# irisデータセット
head(iris)
d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
head(d1)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7

ks.test(x=d1, y="pnorm", mean=mean(d1), sd=sd(d1)) # Kolmogorov-Smirnovの正規性の検定
# Warning message in ks.test(x = d1, y = "pnorm", mean = mean(d1), sd = sd(d1)):
# “ties should not be present for the Kolmogorov-Smirnov test”
#   One-sample Kolmogorov-Smirnov test
#
# data:  d1
# D = 0.1534, p-value = 0.19
# alternative hypothesis: two-sided

シャピロ・ウィルク（Sapiro-Wilk）検定

・サンプルデータから、母集団に正規性があるかどうかを調べたい時に利用する検定
・同じく母集団の正規性を調べる手法として、コルモゴロフ・スミルノフ検定がある

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/Sapiro-WilkTest.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
data = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
display(data.head())
"""
sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)    target
0   5.1 3.5 1.4 0.2 0.0
1   4.9 3.0 1.4 0.2 0.0
2   4.7 3.2 1.3 0.2 0.0
3   4.6 3.1 1.5 0.2 0.0
4   5.0 3.6 1.4 0.2 0.0
0    1.4
1    1.4
2    1.3
3    1.5
4    1.4
Name: petal length (cm), dtype: float64
"""

result = stats.shapiro(data) # Shapiro-Wilkの正規性の検定
print(result) # 統計量、p-value
"""
(0.9549458622932434, 0.05464918911457062)
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/Sapiro-WilkTest.ipynb

# irisデータセット
head(iris)
d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
head(d1)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7

shapiro.test(x=d1) # Sapiro-Wilkの正規性の検定
#
#  Shapiro-Wilk normality test
#
# data:  d1
# W = 0.95498, p-value = 0.05481

フィッシャーの正確確率検定

・独立性の検定
・観測されたデータを水準ごとに分けた時に、その水準間において関係があるかどうかを調べる
・同じく独立性の検定の手法として、カイ二乗検定がある
・カイ二乗検定では、カイ二乗統計量がカイ二乗分布に近似できることから導出されるが、データ数が少ないなどで近似が悪くなる場合には、こちらのフィッシャーの正確確率検定を用いることがある

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/FisherExactTest.ipynb

import pandas as pd
import numpy as np
from scipy import stats

# データ
data = np.array([
    [43, 20],
    [13, 19]
])
"""
        治癒した 治癒しなかった
新薬 43　  20
偽薬 13   19
"""

odds, p=stats.fisher_exact(data) # フィッシャーの正確確率検定
print("オッズ比: ", odds)
print("p-value: ", p)
"""
オッズ比:  3.14230769231
p-value:  0.0147744097701
"""

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/FisherExactTest.ipynb

# データ
data <- matrix(c(43, 20, 13, 19), ncol=2, nrow=2)
#        治癒した 治癒しなかった
# 新薬 43　  20
# 偽薬 13   19

fisher.test(data) # フィッシャーの正確確率検定
#   Fisher's Exact Test for Count Data
#
# data:  data
# p-value = 0.01477
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  1.190240 8.362492
# sample estimates:
# odds ratio
#   3.101933

無相関検定

・2標本データから得られた相関係数を用いて、母相関（真の相関係数）が0であるかそうでないかを調べる手法
・母相関の区間推定については、フィッシャーの $z$ 変換を用いることにより、信頼区間を求めることができる
・無相関検定では、回帰モデルにおける回帰係数の検定の仕組みを用いて、検定統計量を導出できる点を利用している

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/Non-CorrelationTest.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
d1 = data[data["target"] == 0.0] # setosa
print(len(d1))
display(d1.head())
# 50
# sepal length (cm) sepal width (cm)    petal length (cm)   petal width (cm)    target
# 0 5.1 3.5 1.4 0.2 0.0
# 1 4.9 3.0 1.4 0.2 0.0
# 2 4.7 3.2 1.3 0.2 0.0
# 3 4.6 3.1 1.5 0.2 0.0
# 4 5.0 3.6 1.4 0.2 0.0

x = d1["sepal length (cm)"].values
y = d1["sepal width (cm)"].values
r, p = stats.pearsonr(x, y) # 相関係数と無相関検定
print(r) # 相関係数
print(p) # p-value
# 0.746780373264
# 4.75198658015e-10

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/Non-CorrelationTest.ipynb

d1 <- subset(iris, Species=="setosa") # setosa
head(d1)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa

x <- d1$Sepal.Length
y <- d1$Sepal.Width
cor.test(x, y, method="pearson") # 相関係数と無相関検定
#
#   Pearson's product-moment correlation
#
# data:  x and y
# t = 7.6807, df = 48, p-value = 6.71e-10
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  0.5851391 0.8460314
# sample estimates:
#       cor
# 0.7425467

拡張ディッキー-フラー（ADF）検定（Augmented Dickey-Fuller test, ADF test）

・時系列のサンプルデータ（確率過程の生成値）が、単位根過程（非定常過程）であるかどうかを調べたい時に利用する検定
・帰無仮説は「非定常過程である」となり、ADF検定では、非定常過程であることを棄却し、定常過程であるとみなしてよいかどうかを確認する

Python（StatsModels）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/ADFTest_StatsModels.ipynb

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pylab as plt
matplotlib.style.use("ggplot")
from pandas_datareader import data, wb
from statsmodels.tsa import stattools

# データ
start = "2016-08-01"
end = "2017-08-01"
price = data.DataReader("NIKKEI225", "fred", start, end) # 日経平均株価 from セントルイス連邦準備銀行
display(price.head(5))
"""
    NIKKEI225
DATE
2016-08-01  16635.77
2016-08-02  16391.45
2016-08-03  16083.11
2016-08-04  16254.89
2016-08-05  16254.45
"""

# 原系列に対する拡張ディッキー-フラー（ADF）検定
price_tmp = price.dropna()
res_ctt = stattools.adfuller(price_tmp["NIKKEI225"], regression="ctt") # トレンド項あり（２次）、定数項あり
res_ct = stattools.adfuller(price_tmp["NIKKEI225"], regression="ct") # トレンド項あり（１次）、定数項あり
res_c = stattools.adfuller(price_tmp["NIKKEI225"], regression="c") # トレンド項なし、定数項あり
res_nc = stattools.adfuller(price_tmp["NIKKEI225"], regression="nc") # トレンド項なし、定数項なし
print(res_ctt)
print(res_ct)
print(res_c)
print(res_nc)
"""
いずれもp-valueが高く、帰無仮説「非定常過程である」を棄却できない
(-2.6708678219785331, 0.47095101565941583, 1, 245, {'5%': -3.8567110211051521, '1%': -4.4190482250252874, '10%': -3.5682896447228618}, 3022.9312013403046)
(-2.0935231326738535, 0.54972174948713115, 1, 245, {'5%': -3.4285636226572258, '1%': -3.996204153626465, '10%': -3.1376703806237196}, 3023.2796069397737)
(-1.4657871744530102, 0.55036298288220797, 1, 245, {'5%': -2.873410402808354, '1%': -3.4573260719088132, '10%': -2.5730959808413161}, 3023.3326170286091)
(1.3888258295690681, 0.958418587844728, 1, 245, {'5%': -1.9421502633766543, '1%': -2.5749261391087046, '10%': -1.6157794081377657}, 3024.1825324119814)
"""

# 差分系列に対するADF検定
price_tmp = price.dropna()
price_tmp = price_tmp.diff().dropna()
res_ctt = stattools.adfuller(price_tmp["NIKKEI225"], regression="ctt") # トレンド項あり（２次）、定数項あり
res_ct = stattools.adfuller(price_tmp["NIKKEI225"], regression="ct") # トレンド項あり（１次）、定数項あり
res_c = stattools.adfuller(price_tmp["NIKKEI225"], regression="c") # トレンド項なし、定数項あり
res_nc = stattools.adfuller(price_tmp["NIKKEI225"], regression="nc") # トレンド項なし、定数項なし
print(res_ctt)
print(res_ct)
print(res_c)
print(res_nc)
"""
(-18.000023264086924, 0.0, 0, 245, {'5%': -3.8567110211051521, '1%': -4.4190482250252874, '10%': -3.5682896447228618}, 3015.0866820652368)
(-18.033126128850842, 0.0, 0, 245, {'5%': -3.4285636226572258, '1%': -3.996204153626465, '10%': -3.1376703806237196}, 3013.0968669700524)
(-18.045279151781152, 2.6499198050041372e-30, 0, 245, {'5%': -2.873410402808354, '1%': -3.4573260719088132, '10%': -2.5730959808413161}, 3011.6923210889313)
(-17.939018358965527, 5.4207442594331384e-29, 0, 245, {'5%': -1.9421502633766543, '1%': -2.5749261391087046, '10%': -1.6157794081377657}, 3011.9270493988874)
"""

母平均の区間推定

・サンプルデータから、母集団の真の平均が含まれると思われる区間を推定する手法
・信頼係数は、実際に母集団の真の平均が、区間に含まれる確率を表す

Python（Scipy）

GitHub: https://github.com/Gin04gh/samples_py/blob/master/T-StatsConfidenceInterval.ipynb

import pandas as pd
import numpy as np
from sklearn import datasets
from scipy import stats

# irisデータセット
iris = datasets.load_iris()
data = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
display(data.head())
d1 = data[data["target"] == 0.0]["petal length (cm)"] # setosaのPetal.Length
display(d1.head())
# sepal length (cm) sepal width (cm)    petal length (cm)   petal width (cm)    target
# 0 5.1 3.5 1.4 0.2 0.0
# 1 4.9 3.0 1.4 0.2 0.0
# 2 4.7 3.2 1.3 0.2 0.0
# 3 4.6 3.1 1.5 0.2 0.0
# 4 5.0 3.6 1.4 0.2 0.0
# 0    1.4
# 1    1.4
# 2    1.3
# 3    1.5
# 4    1.4
# Name: petal length (cm), dtype: float64

# t-分布による母平均の区間推定
alpha = 0.95
ci = stats.t.interval(alpha, len(d1)-1, loc=np.mean(d1), scale=stats.sem(d1)) # 信頼区間
print(ci)
# (1.4146886741595255, 1.5133113258404745)

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/T-StatsConfidenceInterval.ipynb

d1 <- subset(iris, Species=="setosa")$Petal.Length # setosaのPetal.Length
head(d1)
# 1.4
# 1.4
# 1.3
# 1.5
# 1.4
# 1.7

t.test(d1, mu=0) # t-分布による母平均の検定（帰無仮説「平均は０」）と区間推定
#   One Sample t-test
#
# data:  d1
# t = 59.528, df = 49, p-value < 2.2e-16
# alternative hypothesis: true mean is not equal to 0
# 95 percent confidence interval:
#  1.412645 1.511355
# sample estimates:
# mean of x
#     1.462

母相関の区間推定

・2標本データから、母相関（真の相関係数）が含まれると思われる区間を推定する手法
・信頼係数は、実際に母相関が、区間に含まれる確率を表す
・母相関が無相関がどうかを調べる検定もある

R

GitHub: https://github.com/Gin04gh/samples_r/blob/master/Non-CorrelationTest.ipynb

d1 <- subset(iris, Species=="setosa") # setosa
head(d1)
# Sepal.Length  Sepal.Width Petal.Length    Petal.Width Species
# 5.1   3.5 1.4 0.2 setosa
# 4.9   3.0 1.4 0.2 setosa
# 4.7   3.2 1.3 0.2 setosa
# 4.6   3.1 1.5 0.2 setosa
# 5.0   3.6 1.4 0.2 setosa
# 5.4   3.9 1.7 0.4 setosa

x <- d1$Sepal.Length
y <- d1$Sepal.Width
cor.test(x, y, method="pearson") # 相関係数と無相関検定
#
#   Pearson's product-moment correlation
#
# data:  x and y
# t = 7.6807, df = 48, p-value = 6.71e-10
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  0.5851391 0.8460314
# sample estimates:
#       cor
# 0.7425467