【機械学習_Python】Ridge回帰、Lasso回帰、Elastic Net回帰

Ridge回帰、Lasso回帰、Elastic Net回帰についてのメモ。
本記事は「機械学習」における回帰分析という視点が強め。

Ridge回帰等と区別するため、最小二乗法による一般的な回帰分析を最小二乗回帰分析（ordinary least square：OLS）と呼ぶ。

本記事では Python の sklearn を用い、損失関数の計算式は scikit-learn 0.24.2 の仕様のものを記載している。

【目次】

はじめに① ～最小二乗回帰分析（ordinary least square：OLS）のおさらい～

---

最小二乗回帰分析については次の記事も参照。統計学寄りで書いた記事。
・【統計】単回帰分析 - こちにぃるの日記
・【統計】重回帰分析 - こちにぃるの日記

回帰モデル、予測式及び回帰係数（重み）は次のとおり（上の記事より抜粋）。

これを行列表現すると回帰係数（重み）は次のように求められる。

正規方程式

　　※ ^ t = 転置、^ -1 = 逆行列

Pythonコード

#-- インポート
from sklearn.linear_model import LinearRegression
from sklearn import metrics
import numpy as np
  
#-- 訓練・テスト
mod1 = LinearRegression()
mod1.fit(X_train, y_train)  #-- 訓練データで学習
y_pred1 = mod1.predict(X_test)  #-- テストデータで予測
  
#-- モデル評価
R2sq = mod1.score(X_test, y_test)  #-- 決定係数（寄与率）の算出
MAE  = metrics.mean_absolute_error(y_test, y_pred1)  #--  MAEの算出
RMSE = np.sqrt(metrics.mean_squared_error(y_test, y_pred1)) #-- RMSEの算出
  
#-- 切片、回帰係数（重み）の取得
print("切片 = ", mod1.intercept_)
print("回帰係数 = ", mod1.coef_)

はじめに② ～正則化法の導入～

---

最小二乗回帰分析では、損失関数 を実測値と予測値の差（残差 ）とし、当該関数の最小化を目的とした。

最小二乗回帰分析は機械学習の分野でも「教師あり学習」として広く用いられる手法であるが、複数の特徴量（説明変数）を扱うモデルにおいては次の問題に注意が必要となる。特に「過学習」は将来の予測に重きを置く機械学習の分野で良く耳にするワードである。

・多重共線性（Multicollinearity）
・過学習（Over fitting）

これら問題への対策の１つとして正則化法がある。
式としては、最小二乗回帰分析の損失関数 に正則化項 （ペナルティ項、罰則項）を加えたものを目的関数とする。

正則化項のラムダ はハイパーパラメータであり、ラムダ を０にすれば通常の最小二乗回帰分析となる。
また、 は回帰係数（重み）を用いたLpノルムであり、L1またはL2ノルムが用いられる。損失関数（残差）を小さくするために特徴量（説明変数）を多くすれば、Lpノルムが大きくなり罰則として働く。

L1ノルム（マンハッタン距離）

L2ノルム（ユーグリッド距離）

回帰係数２つのL1ノルムとL2ノルムの2乗は次のとおりである。

なお、Lp正則化を行う場合、特徴量（説明変数）を事前に標準化しておく必要がある。

Pythonコード

from sklearn.preprocessing import StandardScaler
 
SS = StandardScaler()
X_train_SS = SS.fit_transform(X_train)
X_test_SS  = SS.transform(X_test)

Ridge回帰

---

Ridge回帰はL2正則化を用いた方法である。L2ノルム（ユーグリッド距離）の２乗を用いる。

Ridge回帰の特徴は、多重共線性の問題を回避できるところである。
最小二乗回帰分析のセクションで正規方程式を示したが、最小二乗回帰分析では逆行列 が存在しない場合、多重共線性の問題が生じる。

Ridge回帰では、正規方程式に が加わることで逆行列が計算できるようになり、多重共線性の問題が解消される。なお、 は単位行列、 はスカラーでハイパーパラメータである。

Ridge回帰では特徴量選択（変数選択）こそ出来ないが、相関の高いデータセットに対して安定性があり、汎化性の向上が見込める。

Pythonコード

from sklearn.linear_model import Ridge
mod1 = Ridge(alpha=1)
mod1.fit(X_train_SS, y_train_SS)
 
#-- 以降の使い方はLinearRegressionと同じなので割愛

※引数 alpha が本記事内のハイパーパラメータ ：λ（ラムダ）にあたる。コード例の設定値はデフォルト。デフォルトのままで良いなら記載不要。

Lasso回帰（Least absolute shrinkage and selection operator regression）

---

Lasso回帰はL1正則化を用いた方法である。L1ノルム（マンハッタン距離）を用いる。

Lasso回帰の特徴は、いくつかの回帰係数（重み）が０になることであり、この性質から特徴量選択（変数選択）を自動的に行うことができる。スパース推定と言われる。

Lasso回帰による特徴量選択（変数選択）は非常に有用であるが、使用の際は次の性質に注意されたい。

• サンプルサイズn以上の特徴量を選択できない。
• 相関の高い特徴量同士では片方が選択されない可能性がある。重要な特徴量同士の相関が高い場合、特に注意。
• 特徴量間の相関が高いデータセットでは、少しの変動で選択される特徴量が変わる不安定さがあるようである。

回帰係数（重み）の計算方法は現時点で理解できていないため割愛する。
PythonのskleanやRのglmnetではCoordinate Descent法（CD法、座標降下法）が用いられているようである。

次の図は、正則化の概念を示した良く見かける図であり、損失関数（最小二乗誤差）の最適値と制約条件内（水色の円形と四角形）での最適値を示したものである。



制約条件内の最適値が線上にあれば回帰係数（重み）が０になるため特徴量選択ができる。
この制約条件（水色の円形と四角形）の大きさはハイパーパラメータ：λ（ラムダ）で調整できる。

Pythonコード

from sklearn.linear_model import Lasso
mod1 = Lasso(alpha=1)
mod1.fit(X_train_SS, y_train_SS)
  
#-- 以降の使い方はLinearRegressionと同じなので割愛

※引数 alpha が本記事内のハイパーパラメータ ：λ（ラムダ）にあたる。コード例の設定値はデフォルト。デフォルトのままで良いなら記載不要。

Elastic Net回帰

---

Elastic Net回帰はL1正則化とL2正則化の両方を用いた方法である。

L1正則化の「特徴量選択」とL2正則化の「特徴量間の相関の高いデータセットに対する安定性」を組み合わせた方法であり、L1正則化とL2正則化の比率（L1-Ratio）を指定することでモデルの調整ができる。L1-Ratioはハイパーパラメータである。

Ridge回帰やLasso回帰に比べて、Elastic Net回帰が必ずしも良いという訳ではなく、最終的には決定係数やRMSE等の評価指標を用いて最良モデルを選択する。

Elastic Net回帰は式の通り、L1-Ratio=1でLasso回帰と同等となるが、L1-Ratio=0の時はRidge回帰とイコールにならない。
※ただし、Rのglmnetライブラリでは、Elastic Netの式をベースとしてL1-Ratio=1をLasso回帰、L1-Ratio=0をRidge回帰としているようである。
https://cran.r-project.org/web/packages/glmnet/glmnet.pdf
https://cran.r-project.org/web/packages/glmnet/vignettes/glmnet.pdf

Pythonコード

from sklearn.linear_model import ElasticNet
mod1 = ElasticNet(alpha=1, l1_ratio=0.5)
mod1.fit(X_train_SS, y_train_SS)
  
#-- 以降の使い方はLinearRegressionと同じなので割愛

※引数 alphaとl1_ratio が本記事内のハイパーパラメータ ：λ（ラムダ）とρ（ロー、L1-ratio）にあたる。コード例の設定値はデフォルト。デフォルトのままで良いなら記載不要。

bostonデータで各モデルを実行してみる

---

bostonデータを使って最小二乗回帰、Ridge回帰、Lasso回帰、Elastic Net回帰を実行してみる。

from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn import datasets, model_selection, preprocessing, metrics
import pandas as pd
import numpy as np

#-- ロード
ROW = datasets.load_boston()
ADS = pd.DataFrame(ROW.data, columns=ROW.feature_names).assign(target=np.array(ROW.target))

#-- 分割
X_train, X_test, y_train, y_test = model_selection.train_test_split(
        ADS[ROW.feature_names], ADS["target"], test_size=0.3, random_state=1)

#-- 標準化
SS = preprocessing.StandardScaler()
X_train_SS = SS.fit_transform(X_train)
X_test_SS  = SS.transform(X_test)

#-- 回帰分析
lm1 = {"最小二乗回帰":LinearRegression()
       ,"Lasso回帰":Lasso(alpha=1)
       ,"Ridge回帰":Ridge(alpha=1)
       ,"ElasticNet回帰 L1-Ratio=0"  :ElasticNet(alpha=1, l1_ratio=0)
       ,"ElasticNet回帰 L1-Ratio=0.5":ElasticNet(alpha=1, l1_ratio=0.5)
       ,"ElasticNet回帰 L1-Ratio=1"  :ElasticNet(alpha=1, l1_ratio=1)
       }

tmp1 = pd.DataFrame(["決定係数","MAE","RMSE"], columns=["評価指標"])
tmp2 = pd.DataFrame(ROW.feature_names, columns=["特徴量"])

for modnm, modcd in lm1.items(): 
    mod1 = modcd
    mod1.fit(X_train_SS, y_train)
    y_pred1 = mod1.predict(X_test_SS)
    R2sq = mod1.score(X_test_SS, y_test)
    MAE  = metrics.mean_absolute_error(y_test, y_pred1)
    RMSE = np.sqrt(metrics.mean_squared_error(y_test, y_pred1))
    tmp1[modnm] = [R2sq, MAE, RMSE]
    tmp2[modnm] = mod1.coef_

今回の結果をRMSEで比べてみるとRidge回帰が最も良好だった。

また、回帰係数（重み）については、Lasso回帰とElasticNet回帰（L1-Ratio > 0）で特徴量選択がされている。
他、ElasticNet回帰のセクションで触れた通り、L1-Ratio=1ではLasso回帰と同値となるが、L1-Ratio=0ではRidge回帰と一致しない。

参考

---

1. Supervised learning — scikit-learn 0.24.2 documentation
https://www.jstage.jst.go.jp/article/bjsiam/28/2/28_28/_pdf/-char/ja
Rでスパースモデリング：Elastic Net回帰についてまとめてみる - データサイエンティスト(仮)
リッジ回帰とラッソ回帰の理論と実装を初めから丁寧に - Qiita
線形回帰・Lasso(ラッソ)回帰・Ridge(リッジ)回帰を実装してみよう！｜スタビジ

【R】ggplot：対応のあるデータの推移図 ＋ 2つの座標点の角度を求める

対応のあるデータを線でつないで増減をみるグラフを作りたくなったのでメモ。

サンプルデータ

---

library(tidyverse)

ads <- data.frame(id=(1:5)
                  ,x=c(9,6,3,8,7)
                  ,y=c(5,1,2,4,1)
                  ,z=c(8,7,9,5,2))

ads2 <- data.frame()
for (i in (2:4)) {
  ads2 <- union_all(ads2
                   ,data.frame(id=ads[1]
                               ,param=rep(colnames(ads[i]),nrow(ads))
                               ,aval =c(ads[,i])))
}
ads2$param <- factor(ads2$param, level=c("x","y","z"))
ads2$id <- factor(ads2$id)

作図

---

ggplot(data=ads2, mapping=aes(x=param, y=aval, group=id, color=id)) +
  geom_vline(xintercept="x", lty=1) +
  geom_vline(xintercept="y", lty=1) +
  geom_vline(xintercept="z", lty=1) +
  geom_line(size=1.2) +
  annotate("text", x="x", y=ads2[ads2$param=="x","aval"]
           ,label=ads2[ads2$param=="x","id"], hjust=2, size=8) +
  labs(title="対応のあるデータの変化の推移", x="時期", y="値") +
  scale_x_discrete(labels=c("時期1","時期2","時期3")) +
  theme(plot.title  = element_text(size=20)
        ,axis.title = element_text(size=20)
        ,axis.text  = element_text(size=20)
        ,legend.title = element_text(size=20)
        ,legend.text  = element_text(size=20)
        )

2つの座標点の角度を求める

---

角度により増減の度合いを判断したいとする。
すると、コサインを知りたい。
弧度法による計算は次のとおり。atan2（アークタンジェント2）関数を用いる。度数法への変換もできる。
atan2 - Wikipedia

　弧度法 ： ラジアン = atan2(y, x)
　度数法 ： 度 = ラジアン × 180 ÷ 円周率

> #-- ここでは、y=関心のある値の差、x=1として、yの違いを角度で表す。
> ads$ラジアンxy <- atan2(y = ads$y - ads$x, x = 1)
> ads$度xy <- ads$ラジアンxy * 180 / pi
> 
> ads$ラジアンyz <- atan2(y = ads$z - ads$y, x = 1)
> ads$度yz <- ads$ラジアンyz * 180 / pi
> ads
 
  id x y z ラジアンxy      度xy ラジアンyz     度yz
1  1 9 5 8 -1.3258177 -75.96376  1.2490458 71.56505
2  2 6 1 7 -1.3734008 -78.69007  1.4056476 80.53768
3  3 3 2 9 -0.7853982 -45.00000  1.4288993 81.86990
4  4 8 4 5 -1.3258177 -75.96376  0.7853982 45.00000
5  5 7 1 2 -1.4056476 -80.53768  0.7853982 45.00000

【Python】Pandasのmergeとconcat

Pandas の pd.merge と pd.concat の使い方の備忘録。

merge

---

2つのデータフレームを結合する。横結合（join）に対応。

SQLで使われる内部結合（inner join）、外部結合（outer join）、左結合（left join）、右結合（right join）と同じ機能を一通り扱える。

コード例

pd.merge(df1, df2, how="inner", on="USUBJID")

引数（一部）

引数	内容
how =	"inner" で内部結合。 "outer" で外部結合。 "left" で左結合。 "right" で右結合。
on = left_on = right_on = left_index = right_index =	別表参照。
sort =	ソートの有無。Trueでソートする。デフォルトはFalse。
suffixes =	同一列名があった場合、列名末尾に追加する文字を指定。 デフォルトは ('_x' '_y') 。 末尾を_A, _Bとしたければ suffixes = ('_A' '_B') とする。

結合キーの組み合わせ

引数	内容
on = "カラム名"	同名のカラムで結合。
left_on = "カラム名1" ,right_on = "カラム名2"	カラム名をそれぞれ指定して結合。
left_index = True ,right_index = True	インデックス同士で結合。
left_on = "カラム名1" ,right_index = True	左はカラム名、右はインデックスを指定して結合。
left_index = True ,right_on = "カラム名2"	右はカラム名、左はインデックスを指定して結合。

結合方法の組み合わせの例

ADS11 = pd.DataFrame({"ID1":[1,2,3,4,5], "AVAL":[3,5,3,5,6]})
ADS12 = pd.DataFrame({"ID1":[3,4,5,6,7], "ID2":[5,6,7,8,9], "AVAL":[7,8,5,8,7]})
 
#-- ID1で結合。
pd.merge(ADS11, ADS12, how="inner", on="ID1", suffixes=("_11","_12"))
Out[1]: 
   ID1  AVAL_11  ID2  AVAL_12
0    3        3    5        7
1    4        5    6        8
2    5        6    7        5
 
#-- ADS11はID1、ADS12はID2で結合。
pd.merge(ADS11, ADS12, how="inner", left_on="ID1", right_on="ID2")
Out[2]: 
   ID1_x  AVAL_x  ID1_y  ID2  AVAL_y
0      5       6      3    5       7
 
#-- インデックス同士で結合
pd.merge(ADS11, ADS12, how="inner", left_index=True, right_index=True)
Out[36]: 
   ID1_x  AVAL_x  ID1_y  ID2  AVAL_y
0      1       3      3    5       7
1      2       5      4    6       8
2      3       3      5    7       5
3      4       5      6    8       8
4      5       6      7    9       7
 
#-- ADS11はインデックス、ADS12はID1を使って結合
pd.merge(ADS11, ADS12, how="inner", left_index=True, right_on="ID1")
Out[3]: 
   ID1  ID1_x  AVAL_x  ID1_y  ID2  AVAL_y
0    3      4       5      3    5       7
1    4      5       6      4    6       8
 
#-- ADS11はID1、ADS12はインデックスを使って結合
pd.merge(ADS11, ADS12, how="inner", left_on="ID1", right_index=True)
Out[4]: 
   ID1  ID1_x  AVAL_x  ID1_y  ID2  AVAL_y
0    1      1       3      4    6       8
1    2      2       5      5    7       5
2    3      3       3      6    8       8
3    4      4       5      7    9       7

concat

---

２つ以上のデータフレームを結合する。横結合（join）と縦結合（union）の両機能が扱える。

横結合（join）においては、pd.merge と異なり、３つ以上のデータフレームも一度に扱える特徴を持つ。 一方、結合方式は内部結合（inner join）と外部結合（outer join）の２択になるため、一長一短。使い分けになる。

コード例

pd.concat([df1, df2], axis=0, join="outer", keys=["Group1", "Group2"])
pd.concat([df1, df2, df3], axis=0, join="outer", keys=["Group1", "Group2"])

引数（一部）

引数	内容
axis =	0 で縦結合、1 で横結合。
join =	"outer" で全ての列／行を結合。 "inner" で共通する列／行のみ結合。
ignore_index =	Trueでインデックスを振りなおす。 デフォルトは False で、結合前のインデックスをそのまま使う。
keys =	リスト型でマルチインデックス内のキー名を指定する。 例：["Group1", "Group2"] ※マルチインデックスの設定例を参照。
names =	リスト型でインデックス名を指定する。 例：["index1", "index2"] ※マルチインデックスの設定例を参照。
sort =	ソートの有無。Trueでソートする。デフォルトはFalse。

結合方法の組み合わせの例

ADS11 = pd.DataFrame({"COL1":[1,3], "AVAL":[3,5]}, index=[1,2])
ADS12 = pd.DataFrame({"COL2":[3,5], "AVAL":[7,8]}, index=[3,2])
 
#-- 縦結合。全ての列を結合。
pd.concat([ADS11, ADS12], axis=0, join="outer")
Out[1]: 
   COL1  AVAL  COL2
1   1.0     3   NaN
2   3.0     5   NaN
3   NaN     7   3.0
2   NaN     8   5.0
 
#-- 縦結合。共通列のみ結合。
pd.concat([ADS11, ADS12], axis=0, join="inner")
Out[2]: 
   AVAL
1     3
2     5
3     7
2     8
 
#-- 横結合。全ての列を結合。
pd.concat([ADS11, ADS12], axis=1, join="outer")
Out[3]: 
   COL1  AVAL  COL2  AVAL
1   1.0   3.0   NaN   NaN
2   3.0   5.0   5.0   8.0
3   NaN   NaN   3.0   7.0
 
#-- 横結合。共通列のみ結合。
pd.concat([ADS11, ADS12], axis=1, join="inner")
Out[4]: 
   COL1  AVAL  COL2  AVAL
2     3     5     5     8

マルチインデックスの設定例

#-- 縦結合
pd.concat([ADS11, ADS12]
            ,axis = 0, join = "outer"
            ,keys = ["Group1", "Group2"]
            ,names = ["Index1","Index2"]
            ,sort   = True
    )
 
Out[1]: 
               AVAL  COL1  COL2
Index1 Index2                  
Group1 1          3   1.0   NaN
       2          5   3.0   NaN
Group2 3          7   NaN   3.0
       2          8   NaN   5.0
 
#-- 横結合
pd.concat([ADS11, ADS12]
            ,axis = 1, join = "outer"
            ,keys = ["Group1", "Group2"]
            ,names = ["Index1","Index2"]
            ,sort   = True
    )
 
Out[2]: 
Index1 Group1      Group2     
Index2   COL1 AVAL   COL2 AVAL
1         1.0  3.0    NaN  NaN
2         3.0  5.0    5.0  8.0
3         NaN  NaN    3.0  7.0

【Python】DataFrame型とSeries型

PandasのDataFrame型とSeries型の備忘録

サンプルデータ

---

import pandas as pd
group1 = [4, 9, 10, 11, 12, 13, 13, 15, 18, 18, 20]
group2 = [16, 16, 17, 19, 19, 22, 22, 23, 25]

DataFrame型

---

DataFrame（データフレーム）型は表（テーブル）形式の２次元データ。

addf = pd.DataFrame({"AVAL":group1 + group2
                    ,"TRT01A":["group1"] * len(group1) + ["group2"] * len(group2) })

↓中身

addf.head(5)
Out[1]: 
   AVAL  TRT01A
0     4  group1
1     9  group1
2    10  group1
3    11  group1
4    12  group1
 
type(addf)
Out[2]: pandas.core.frame.DataFrame

Series型

---

Series（シリーズ）型は１列または１行の1次元データ。

adsr = pd.Series(data=group1, index=range(1,len(group1)+1))

↓　中身

adsr.head(5)
Out[1]: 
1     4
2     9
3    10
4    11
5    12
dtype: int64
 
type(adsr2)
Out[2]: pandas.core.series.Series

DataFrame型データを１列取り出すと１次元データになるので、Series型と認識される。

adsr2 = addf.loc[addf["TRT01A"]=="group1", "AVAL"]

↓　中身

adsr2.head(5)
Out[86]: 
0     4
1     9
2    10
3    11
4    12
Name: AVAL, dtype: int64

type(adsr2)
Out[87]: pandas.core.series.Series

【Python】Classの定義（イニシャライザ、継承も含めて）

クラスを試してみる。

クラス、インスタンス、インヘリタンスの定義

---

・クラス（class）：
設計図と言われる。このクラスから実体（インスタンス）を生成する。
クラス自体に実体はない（データを持たない）。

・インスタンス（instance）：
日本語で実体。
クラスにデータを与えて生成されるもの。
実体を生成することをインスタンス化と言う。

・インヘリタンス（inheritance）：
日本語で継承。
親から子へクラスの内容を引き継がせる機能。

Pythonでの実装（クラス定義、インスタンス化）

---

Pythonでのクラス定義は class クラス名 : と書き、インデントを下げて関数を定義していく。
クラス内の関数はメゾッドと呼ばれ、第一引数で自分自身を定義する必要がある。第一引数の名称は self とする慣習がある。

また、次のコード例で示す def __init__(self, val1...) は特別な関数でイニシャライザと呼ばれる（初期値を与える処理）。

#-----------------------------
#-- クラス定義
#-----------------------------
class MyClass1: 
    def __init__(self, val1, val2):  #-- 初期値を与える
        self.val1 = val1
        self.val2 = val2
    
    def 足し算(self):  #-- メゾッド
        x = self.val1 + self.val2
        print ("足し算：{0}+{1}={2}".format(self.val1,self.val2,x))
 
    def 引き算(self):  #-- メゾッド
        x = self.val1 - self.val2
        print ("引き算：{0}-{1}={2}".format(self.val1,self.val2,x))
 
    def 掛け算(self):  #-- メゾッド
        x = self.val1 * self.val2
        print ("掛け算：{0}×{1}={2}".format(self.val1,self.val2,x))

次に、インスタンス名 = クラス名 (引数) とすることで、定義したクラスに値を与えてインスタンス化する。
そして、インスタンス名.メゾッド名() とすることで、クラス内に定義したメゾッドを実行する。

#-----------------------------
# インスタンス化
#-----------------------------
res1 = MyClass1(1, 2)  #-- インスタンス化（MyClass1に値を与えて実体化する）
  
#-----------------------------
# クラス定義内のメゾッドを実行
#-----------------------------
res1.足し算()
res1.引き算()
res1.掛け算()
 
--- out ---
足し算：1+2=3
引き算：1-2=-1
掛け算：1×2=2

上のコード例の通り、クラス定義により様々な機能を１つの括りでまとめることができる。
大規模な開発になるほど可読性やメンテナンス性などが向上する可能性がある。

Pythonでの実装（継承）

---

クラスでは親クラス→子クラスに値を引き継がせる継承機能を扱うことができる。

子クラスは class クラス名 (親クラス名) : とする。
この時、子クラス内で super().親クラスのメゾッド(引数名) とすれば、親クラスのメゾッドを呼び出すことができ、同一処理をわざわざ再定義する必要がない。

せっかくなので、これまでの自身の記事を参考にして、
スチューデントのT検定 と Wilcoxonの順位和検定 を実行できる簡単なパッケージを作ってみる。（※検定名のところに記事のリンク貼っている）

class MyStat1:
    def __init__(self, group1, group2):
        self.group1 = group1
        self.group2 = group2
        self.n1 = len(group1)  #-- グループ1の例数
        self.n2 = len(group2)  #-- グループ2の例数
        self.N  = self.n1 + self.n2  #-- 全体の例数
 
class StudentTtest(MyStat1):
    def __init__(self, group1, group2):
        super().__init__(group1, group2)  #-- 親クラスのイニシャライザを継承
        print("継承しました")  #-- 継承したか見たいだけのプリント
        
    def Diffmean(self):
        import numpy as np
        m1 = np.mean(group1)
        m2 = np.mean(group2)
        return m1 - m2
        
    def TestResult(self):
        import numpy as np
        from scipy import stats
        s1 = np.sum((group1 - np.mean(group1))**2)
        s2 = np.sum((group2 - np.mean(group2))**2)
        Df = self.N - 2
        se = np.sqrt((s1 + s2) / Df)
        T  = self.Diffmean() / (se * np.sqrt(1/self.n1 + 1/self.n2))
        P  = stats.t.cdf(x=T, df=Df)*2
        print(" ----- Student's T-test -----\n"
              ,"T統計量：{0}\n 自由度：{1}\n p値：{2}\n".format(T, Df, P)
              ,"----------------------------")
        return T, Df, P
    
class WilcoxonRankSumTest(MyStat1):
    def __init__(self, group1, group2, ties=True):
        super().__init__(group1, group2)  #-- 親クラスのイニシャライザを継承
        self.ties = ties  #-- Wilcoxonの順位和検定はさらにタイ補正有無の引数を与えてみる
        print("継承しました")  #-- 継承したか見たいだけのプリント
        
    def TestResult(self):
        import numpy as np
        from scipy import stats
        import pandas as pd
        
        #-- データフレーム化
        ads = pd.DataFrame({"AVAL":group1 + group2
                            ,"TRT01AN":["group1"] * self.n1 + ["group2"] * self.n2
                            })
        #-- 全体順位
        ads["Rank"] = ads["AVAL"].rank()
        
        #-- タイ補正（同一順位の数で補正）
        adfrq = ads["AVAL"].value_counts()
        wil_t = np.sum(adfrq**3 - adfrq)
        
        if self.ties == False:
            wil_t = 0  #-- タイ補正しない場合、0に上書き。
        
        #-- 順位和        
        sumrank1 = float(np.sum(ads.loc[ads["TRT01AN"]=="group1", ["Rank"]]))
        sumrank2 = float(np.sum(ads.loc[ads["TRT01AN"]=="group2", ["Rank"]]))
        
        #-- 実測、期待値、分散        
        if sumrank1 < sumrank2:
            wil_W = sumrank1
            wil_E = self.n1 * (self.n1 + self.n2 + 1) / 2
        else:
            wil_W = sumrank2
            wil_E = self.n2 * (self.n1 + self.n2 + 1) / 2
            
        val1 = self.n1 * self.n2 * (self.n1 + self.n2 + 1) / 12
        val2 = self.n1 * self.n2 / (12 * (self.n1 + self.n2) * (self.n1 + self.n2 - 1))
        wil_V = val1 - val2 * wil_t
        wil_Z = (wil_W - wil_E) / np.sqrt(wil_V)
        
        wil_P = stats.norm.cdf(wil_Z)*2
        
        print("----- Wilcoxon Rank Sum Test -----\n"
              ,"順位和：{0} (group1), {1} (group2)\n".format(sumrank1, sumrank2)
              ,"W：{0}\n".format(wil_W)
              ,"E[W]：{0}\n".format(wil_E)
              ,"Var[W]：{0}\n".format(wil_V)
              ,"タイ補正の係数：{0}\n".format(wil_t)
              ,"\n"
              ,"Z統計量：{0}\n".format(wil_Z)
              ,"p値：{0}\n".format(wil_P)
              ,"----------------------------------")
        return wil_Z, wil_P

出来上がったクラスを試してみる。
サンプルデータは Wilcoxonの順位和検定 の記事で扱ったものと同じ。

group1 = [4, 9, 10, 11, 12, 13, 13, 15, 18, 18, 20]
group2 = [16, 16, 17, 19, 19, 22, 22, 23, 25]

T検定（スチューデント）

res1 = StudentTtest(group1, group2)
  
res1.Diffmean()
Out[1]: -6.888888888888889
 
res1.TestResult()
継承しました
 ----- Student's T-test -----
 T統計量：-3.761258193347809
 自由度：18
 p値：0.0014296297831112022
 ----------------------------
Out[1]: (-3.761258193347809, 18, 0.0014296297831112022)

親クラスのN数を self.n1, self.n2, self.N として子クラスで使っているが、問題なく使用できている。
また、StudentTtest クラス内の TestResult メゾッドで、同一クラス内の Diffmean メゾッドを self.Diffmean() で呼び出しているが、こちらも問題はなく動作している。
念のため検算。

from scipy import stats
stats.ttest_ind(group1, group2, equal_var = True)
Out[1]: Ttest_indResult(statistic=-3.7612581933478086, pvalue=0.001429629783111206)

Wilcoxonの順位和検定（タイ補正あり）

res2 = WilcoxonRankSumTest(group1, group2, ties=True)
 
res2.TestResult()
継承しました
----- Wilcoxon Rank Sum Test -----
 順位和：77.0 (group1), 133.0 (group2)
 W：77.0
 E[W]：115.5
 Var[W]：172.59868421052633
 タイ補正の係数：30
 
 Z統計量：-2.930501777999723
 p値：0.00338415057040675
 ----------------------------------
Out[1]: (-2.930501777999723, 0.00338415057040675)

Wilcoxonの順位和検定（タイ補正なし）

res3 = WilcoxonRankSumTest(group1, group2, ties=False)
 
res3.TestResult()
継承しました
----- Wilcoxon Rank Sum Test -----
 順位和：77.0 (group1), 133.0 (group2)
 W：77.0
 E[W]：115.5
 Var[W]：173.25
 タイ補正の係数：0
 
 Z統計量：-2.9249881291307074
 p値：0.0034446936330652616
 ----------------------------------
Out[1]: (-2.9249881291307074, 0.0034446936330652616)

念のため検算。余談だが、stats.ranksumsはタイ補正なしのようだった。

from scipy import stats
stats.ranksums(group1, group2)
Out[2]: RanksumsResult(statistic=-2.9249881291307074, pvalue=0.0034446936330652616)

本記事はここまで。

【Python】関数の作成（def）

関数は「 def 関数名 (引数名) 」とし、インデントを下げてコードを書く。
戻り値は「return 変数名」と書く。なくても良い。

def def_test1(i, j):
    x = i + j
    return x

実行は関数名（引数名）と書く。

val1 = def_test1(1, 2)
print(val1)
  
-- out --
3

---

もう少し規模大きくして、t 検定（Student）の t 統計量を計算する関数を作ってみる。
計算式は 過去記事 参照。

#-- 関数：T統計量を算出してみる
def ttest1(group1, group2):
    import numpy as np
    m1 = np.mean(group1)
    m2 = np.mean(group2)
    n1 = len(group1)
    n2 = len(group2)
    N  = n1 + n2
    s1 = np.sum((group1 - np.mean(group1))**2)
    s2 = np.sum((group2 - np.mean(group2))**2)
    se = np.sqrt((s1 + s2) / (N - 2))
    t  = (m1 - m2) / (se * np.sqrt(1/n1 + 1/n2))
    return t

実行

group1 = [8,6,7,6,8,9]
group2 = [3,4,3,5,4,6]

ttest1(group1, group2)
 
-- out --
4.608176875690326

stats.ttest_ind関数で答え合わせ。

from scipy import stats
stats.ttest_ind(group1, group2, equal_var = True)
   #-- equal_var=True: Student, False: Welch
 
-- out --
Ttest_indResult(statistic=4.608176875690326, pvalue=0.0009679292740393189)

【Python】条件分岐（if文）やループ処理（while文、for文）

Pythonの基礎。
条件分岐やループ処理について。

if 文

---

if 文は他の言語と大体同じ。ただし、else if は elif と書く。

a=1
if a==0:
    print(0)
elif a==1:
    print(1)
else:
    print(2)

while 文

---

while 文は特記なし。よくある形。

cnt=3
while cnt != 0:
    print(cnt)
    cnt -= 1
    if cnt==0:
        print("End!!")
 
-- out --
3
2
1
End!!

continue で強制的に次のループに回す、break で強制的にループを終了できる。

loop_lst=["a", "b", "c", "d", "e"]
cnt = -1
 
while cnt <= len(loop_lst):
    cnt += 1
    print(cnt)
    
    if cnt == 0:
        print("continue!!")
        continue  #-- 以降処理やらず、次のループ処理に回る。
    
    if cnt == 3:
        print("break!!")
        break  #-- ループを終了する
 
-- out --
0
continue!!
1
2
3
break!!

for 文

---

for の対象の与え方はいくつかある。

・レンジを与える方法

for i in range(0, 10, 2):
    print(i)
 
-- out --
0
2
4
6
8

・リストを与える方法

loop_lst=["a","b","c","d","e"]
 
for i in loop_lst:
    if i == "a":
        print("continue!!")
        continue  #-- 以降処理やらず、次のループ処理に回る。
    
    if i == "d":
        print("break!!")
        break  #-- ループを終了する
        
    print(i)
  
-- out --
continue!!
b
c
break!!

・リストを与える方法（インデックスと値を取得）

for_lst1 = ["a", "b", "c", "d", "e"]
 
for index1, val1 in enumerate(for_lst1): #-- enumerateで囲う
    print("{0} {1} {2}".format(index1, val1, 100 + index1))
 
-- out --
0 a 100
1 b 101
2 c 102
3 d 103
4 e 104

・リストを与える方法（入れ子のリスト）

for_lst2_1 = ["a", "b", "c"]
for_lst2_2 = ["あ", "い", "う"]
 
for i, j, k in [for_lst2_1, for_lst2_2]:
    print("要素１：{0}、要素２：{1}、要素３：{2}".format(i, j, k))
 
-- out --
要素１：a、要素２：b、要素３：c
要素１：あ、要素２：い、要素３：う
 
for i, j in zip(for_lst2_1, for_lst2_2):
    print("リスト１：{0}、リスト２：{1}".format(i, j))
 
-- out --
リスト１：a、リスト２：あ
リスト１：b、リスト２：い
リスト１：c、リスト２：う

・辞書を与える方法（キーを取得）

for_dict1 = {"key1":"a", "key2":"b", "key3":"c"}

#-- keysメゾッドでキーを取得
for key in for_dict1.keys():
    print("キー：{0}".format(key))
 
-- out --
キー：key1
キー：key2
キー：key3

・辞書を与える方法（値を取得）

for_dict1 = {"key1":"a", "key2":"b", "key3":"c"}
  
#-- valuesメゾッドで値を取得
for val in for_dict1.values():
    print("値：{0}".format(val))
 
-- out --
値：a
値：b
値：c

・辞書を与える方法（キーと値を取得）

for_dict1 = {"key1":"a", "key2":"b", "key3":"c"}
 
#-- itemsメゾッドでキーと値を取得
for key, val in for_dict1.items():
    print("キー：{0}、値：{1}".format(key, val))
 
-- out --
キー：key1、値：a
キー：key2、値：b
キー：key3、値：c