2019-09-29

統計検定2級(CBT)合格記

統計学統計検定

こんにちは!

company-rapperです!

先日、統計検定2級(CBT)を受験し、合格しましたので、報告したいと思います。

統計検定2級とは?
どうして受けたのか?
受験までに勉強した内容・教材
- 受験2ヶ月前〜1ヶ月前
- 受験3週間前頃〜
受験当日
ご参考: 過去問題のcompany-rapperの正答数と本番試験結果
おわりに

統計検定2級とは?

統計検定は、文字通り、統計学に関する知識やそれらをもとに問題を解くことができるかをテストする検定試験です。

www.toukei-kentei.jp

2級は大学基礎課程のレベルですので、大学1,2年生の時に開講される統計学の授業のレベルとざっくりと理解できます。詳細は上のリンクでご確認ください。

さて、統計検定は年に2回一般会場での試験があります。それに加えて、CBT(Computer Based Testing)と呼ばれるPCを使った試験形式もあります。

こちらは、テストセンターや提携しているパソコン教室などにおいてパソコンを使って検定試験を受けるタイプのもので、(日程は検定会場によりマチマチではあるものの、試験会場を選ばなければ)ほぼ毎日受験することができます。

どうして受けたのか?

データ分析をし始めて「このデータ分析をしている時に感じるモヤモヤは統計学がわかっていないからだ!!」と気づきました。

実は、私、大学の統計学の授業は受けませんでした。なんとなく、最初の授業を受けた時にめんどくさい感じがしたのです(あの時受けていたらなー)

統計学を基本から学び直そう!と土日を使って勉強をし始めたのですが

勉強に時間がかかる!!!!
どこまで勉強すればいいのかがいまいちわからない!!!!

という不満が出てきました。

まず、1番目の「時間がかかる」統計学は未知の分野だったのと統計学独特の考え方に慣れるまで時間がかかりました。せっかく時間をかけたのだから何かの証が欲しい・・・

2番目の「勉強の範囲」ですが、何がわかれば/何を解ければ一通り基礎はできたと言えるのかがいまいち不透明でした。

この両方の不満を解消するのが統計検定だったのです。統計検定が取れれば、なにより、統計学の基礎を勉強したという証がもらえます。プラスして、統計検定に出てくる問題が解けるという一つのクライテリアが生まれるので、自分がちゃんとわかっているかわかっていないかの判断が付きやすくなりました。

さらに、検定の受験日を決めることで、その日に向かって準備しないといけないので、良い勉強のペースメーカーとなりました。

受験までに勉強した内容・教材

受験2ヶ月前〜1ヶ月前

「赤本」で勉強しました

統計学入門 (基礎統計学?)

作者: 東京大学教養学部統計学教室
出版社/メーカー: 東京大学出版会
発売日: 1991/07/09
メディア: 単行本
購入: 158人クリック: 3,604回
この商品を含むブログ (79件) を見る

第6章以降は練習問題も真面目に解きました。割と骨のある問題が多かったかと思います。

受験3週間前頃〜

過去問を解きました。過去問は公式問題集を使いました

日本統計学会公式認定統計検定 2級公式問題集[2016〜2018年]

作者: 日本統計学会
出版社/メーカー: 実務教育出版
発売日: 2019/03/20
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

この1冊で6回分収録されています。試験までに一通り解くのと、解説を読んで知識が怪しいところや間違えそうなところはノートに書いていくことをしました。

f:id:company-rapper:20190928170458j:plain

試験の申し込みは受験の2週間前くらいにしました。過去問を解いて、この正答率なら受かるんじゃないかと思えたタイミングです。CBTだと、自分の受かりそうと思えるベストタイミングで試験日程が組めるので良いと思います。

過去問を解くにつれて、「赤本」には詳しく書いていないが、統計検定2級の範囲に入っている(問題として登場する)項目があることがわかってきました(下記)

ラスパイレス価格指数
データ抽出の方法
ローレンツ曲線・ジニ係数
フィッシャーの3原則
コレログラム
一元配置分散分析(ANOVA)

1〜5に関しては、統計検定2級の公式問題集の解説で勉強しました。

6の一元配置分散分析については

こちらのpdf(読めば必ずわかる分散分析の基礎)がわかりやすかったです。

受験当日

ついに受験当日です。お手製ノートを見直して、見直して、試験に臨みました。いくつか、失敗談がございますので、書きます。

注意点1: 電卓

統計検定では電卓が必須となっています。

f:id:company-rapper:20190929210237j:plain

私、とりあえず、イトーヨーカドーで簡単な電卓を買いました。その電卓を使っていざ問題演習・・・と思いきや、、、、

ルートがない!!!!!

なるべく安いものをと思い過ぎたせいか、ルート機能が付いていない電卓を買ってしまいました。

仕方ないので、

シャープ手帳タイプ電卓高級感を演出するアクリルパネル採用。ブラック系 EL-WA21X

出版社/メーカー: シャープ(SHARP)
発売日: 2015/07/10
メディア: オフィス用品
この商品を含むブログを見る

こちらをAmazonで買いました。特に、テスト前の持ち物チェックで何も言われなかったので大丈夫な電卓だと思います。

注意点2: テストセンター会場

私の受験会場は

otc.odyssey-com.co.jp

こちらでした。

横浜そごうの隣のビルの26Fです。こちらの会場、ビルの一室で、なんと、試験開始(確か)15分前にならないと入室できないというところでした。

26Fのエレベーターホールには特に椅子もなく、1Fのカフェも満員で、仕方なくこのビルの最上階のレストラン街の椅子で試験直前の最終チェックをしました。

が、ママ会が終わってちょうどレストランから出てきた親子連れが同じフロアでお話ししており、あまり集中できず。休みの日は喫茶店も混んでいるのでなかなか集中できる場所を確保するのは大変とわかりました。

注意点3: CBT試験特有の問題構成

統計検定2級では、ペーパーテストですと、大問ひとつに対して問題が2〜3個あるというのが一般的です。この大問が18個程度あるかと思います。

しかしながら、CBTテストは、問題データベースからランダムに問題を抽出する関係上、大問一つに対して小問が1つしかありません!!

これはかなり痛手でした。全問題数はペーパーテストとほぼ同じと思いますので、一問一問大問を理解しないといけないため、過去問を解いている時よりも解くスピードが落ちます。問題を解くごとに問題設定が変わるので頭もいつも以上に疲れます。

時間配分のペースがよくわからなくなり、何問かよくわからない問題は適当にマークし、時間終了となりました。

ペーパーテストですと合格基準が70点となっている一方で、CBTは60点と低めになっているのはこの理由(実質的に難易度が上がるから)かと思います。

ご参考: 過去問題のcompany-rapperの正答数と本番試験結果

ご参考までに、過去問チャレンジ結果と本番の時の結果です。過去問題集には点数配分は書いていないので、正答率でおよその点数を把握しました。

2016/6 31/34問(正答率91%)
2016/11 27/34問(正答率79%)
2017/6 33/35問(正答率94%)
2017/11 28/34問(正答率82%)
2018/6 23/34問(正答率68%)
2018/11 27/34問(正答率79%)

本番: 74点/100点 (合格点60点)

セクション分析(正答率)

1変数・2変数記述統計分野 90%
データ収集・確率・分布の分野 91%
推定・検定・線形モデルの分野 46% (くやしい。もっと取りたかった)

おわりに

昨今のデータ分析に対する盛り上がりに伴い、統計学の重要性はますます高まるばかりです。しかしながら、高校数学では統計は必須ではないですし、私自身、大学の統計学は受けませんでした。

私のように、理系なんだけれども、実は統計学はちゃんと勉強してなかったんだよね・・・でも学び直したい!という方々にはオススメの資格だと思います。

以上、company-rapperでした!

2019-09-21

インターネット速度比較 : GIGA PRIZE(積和不動産シャーメゾンインターネット) vs ソネット(So-net 光プラスマンション IPv6プラスオプション)

ネット

こんにちは

company-rapperです!

先日、アパート各戸に光回線が通りました。いわゆる「無料インターネット付き部屋」になりました。

今住んでいるアパートで空き部屋が埋まらないので、インターネット付きにして価値を高めたいという大家さんの考えがあるのかもしれません。

こんな感じで

f:id:company-rapper:20190921223040j:plain — LANコンセント

電話回線のコンセント部分に有線LANソケットを備えた無線LANの小さな機械が取り付けられます。(光コンセントはもともと付いていました。)

私の部屋まではLANケーブルを壁伝いに引いていたので、おそらく、アパート全体で1回線契約していて、そこからハブで分岐し、各戸までLANで繋いでいるものと思います。工事中、電源を別途繋げている様子はなかったので、Power over Ethernetを使ってLANケーブルで給電しているのだと思います。

コンセントには有線LAN端子が一つあり、無線LANは2.4GHz帯(802.11b/g/n)、5GHz帯(802.11ac)の両方対応です。

工事自体は40分程度で終わりました。

さて、すでに私の部屋にはBフレッツ(So-net)が開通しているので、どちらが早いかを試してみることにしました!!

土曜日の夜9時にテストしました。どちらもWi-Fi(5GHz帯, 802.11ac)を使用しています。

Googleの速度テスト結果です

f:id:company-rapper:20190921223702p:plain — So-net 光プラスマンション (IPv6オプション)

f:id:company-rapper:20190921223749p:plain — GIGA PRIZE(無料インターネット)

次に、BNRスピードテストの結果です

f:id:company-rapper:20190921223926p:plain — So-net 光プラスマンション (IPv6オプション)

f:id:company-rapper:20190921224054p:plain — GIGA PRIZE(無料インターネット)

最後に、FAST (Netflix)での結果です

f:id:company-rapper:20190921224350p:plain — So-net 光プラスマンション(IPv6オプション)

f:id:company-rapper:20190921224439p:plain — GIGA PRIZE(無料インターネット)

さて、まとめます

ダウンロード速度(単位:Mbps)
種類	Google	BNR	FAST
So-net	123.6	41.98	70
無料インターネット	32.6	34.83	170

結局、速度測定サイトによりバラツキが大きくどっちが早いのかよくわかりませんでした。

無料インターネットの方は、前述の通り、アパート全体で1回線分を共有しているので、他の入居者の方々の使い具合で速度はかなり変わるものと思います。

そして、無料インターネットの方は、特に説明書などは渡されませんでした。このため、無線LANの親機側の設定は(おそらく)できないものと思います。あくまでもライトユーザー向けと感じます。

とりあえず、So-netは解約はしないことにして(解約金が結構かかるのです・・・)、無料インターネットはテレビにつなげることとしました。最近の機種はテレビでYouTubeが見れていいですね。

では、company-rapperでした!

2019-09-14

sklearn.preprocessing.PolynomialFeatures のカラム名を取得する〜Boston Housingを例にして〜

機械学習 python scikitlearn

こんにちは、company-rapperです!

前回のポスト

company-rapper.hatenablog.jp

で、sklearn.preprocessing.PolynomialFeatures()の仕様を調査しました。

今回は、その応用編として、PolynomialFeatures()で多項式化させたデータセットを回帰分析して、回帰係数がどのカラム名に対応しているかを出力してみたいと思います。

コード(Jupyter notebook形式)はGithubにアップしています。

github.com

今回は、例としてBoston Housingのデータセットを使用します。これは1970年の国勢調査で得られたボストンの住宅価格に関するデータです。scikit-learnのサンプルデータセットに収録されています。

プログラムの流れとしては

ライブラリインポート
データインポート
今回はBoston Housingを使用
グリッドサーチ
グリッドサーチ内では
StandardScaler()で標準化する
PolynomialFeatures()で多項式カラムを作る
ElasticNet()で回帰させる
を行なっています。
ベストモデルで係数出力(係数が大きい順に出力)
この例では、すなわち、ElasticNetで得られた回帰係数がどの特徴量に対応しているかをランキング形式で出力します。

です

まずはライブラリをインポートします。

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import GridSearchCV
from sklearn.datasets import load_boston

次に、データセット *1をインポートし、学習データとテストデータに分けます

boston = load_boston()
# sklearn.utils.BunchからDataFrame形式へ変換
df_prices = pd.DataFrame(boston.target, columns=["MEDV"])
df_features = pd.DataFrame(boston.data, columns=boston.feature_names)

df_data = pd.concat([df_prices,df_features],axis=1)

# 価格が50以上の場合は50になってしまっているようなので、データを削除
# (「Boston Housing Prices - Evaluation & Validation」https://www.kaggle.com/samratp/boston-housing-prices-evaluation-validation)
df_data = df_data[df_data["MEDV"]<50]

df_features_ = df_data.drop(["MEDV"],axis=1)
df_prices = df_data["MEDV"]
X_train, X_test, y_train, y_test = train_test_split(df_features_, df_prices, test_size = 0.2, random_state=5)

次に、グリッドサーチ*2でハイパーパラメータを探索します

alpha_params = np.logspace(-3, 3, 7)
l1_ratio_params = np.arange(0.1, 1.0, 0.1)
 
steps = [
    ('ss', StandardScaler()),
    ('pf', PolynomialFeatures()),
    ('en', ElasticNet())
]
pipeline = Pipeline(steps=steps)
    
paramters = {
    'pf__degree': [1,2,3],
    'en__alpha': alpha_params,
    'en__l1_ratio': l1_ratio_params,
    'en__fit_intercept': [False]
}
 
model_CV = GridSearchCV(
    estimator = pipeline,  
    param_grid = paramters,
    cv = 5,
    n_jobs = -1, 
)
 
model_CV.fit(X_train, y_train)

best_param = model_CV.best_params_
best_estimator = model_CV.best_estimator_
 
pred_train = best_estimator.predict(X_train)
r2score_train = r2_score(y_train, pred_train)
rmse_train = np.sqrt(mean_squared_error(y_train, pred_train))

pred_test = best_estimator.predict(X_test)
r2score_test = r2_score(y_test, pred_test)
rmse_test = np.sqrt(mean_squared_error(y_test, pred_test))

print("結果")
print('RMSE (train) : {:.3f}'.format(rmse_train))
print('RMSE (test) : {:.3f}'.format(rmse_test))
print('R2 score (train): %.3f' % r2score_train)
print('R2 score (test) : %.3f' % r2score_test)
print('clf.best_params_', best_param)

出力は下記のようになります

RMSE (train) : 2.585
RMSE (test) : 2.798
R2 score (train): 0.889
R2 score (test) : 0.885
clf.best_params_ {'en__alpha': 0.1, 'en__fit_intercept': False, 'en__l1_ratio': 0.9, 'pf__degree': 2}

この出力から、ベストなモデルとなるハイパーパラメータはElasticNetのα=0.1, l1_ratio=0.9, PolynomialFeaturesの次元が2の時と得られました。R2 scoreで見ると、trainが0.889, testが0.885と過学習している心配はなさそうです。(※ElasticNetのfit_intercept=Falseとしているのは、PolynomialFeaturesですでにバイアス項を生成させているためです。)

さて、ここからが本番です。

ElasticNetは回帰型の分析ですから、回帰式の係数が分析の結果として得られています。この例では、

best_estimator.named_steps["en"].coef_

でその係数を得られますが、これはnumpy.ndarray形式であり、出力すると

array([ 1.75804749e+01, -2.04597703e+00, -4.86938426e-01, -9.23279354e-01, -0.00000000e+00, -4.54373705e-01, 2.50403603e+00, -7.75038911e-01, -7.67195961e-01, -0.00000000e+00, -9.50376753e-01, -5.91575503e-01, 5.94885088e-01, -2.48676454e+00, 1.18679136e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -0.00000000e+00, 1.38849955e-02, -0.00000000e+00, 0.00000000e+00, -0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -1.02503909e-01, -0.00000000e+00, 1.54255789e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.10347954e-01, -0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.46273726e-01, -0.00000000e+00, -0.00000000e+00, 5.59900185e-01, -0.00000000e+00, 0.00000000e+00, -0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -0.00000000e+00, -0.00000000e+00, -0.00000000e+00, 1.70831778e-01, -2.04649122e-01, -4.64448865e-01, -0.00000000e+00, 4.09789119e-01, -0.00000000e+00, -0.00000000e+00, -0.00000000e+00, -0.00000000e+00, 0.00000000e+00, -0.00000000e+00, -4.38022921e-01, -0.00000000e+00, 4.22304702e-01, -0.00000000e+00, -0.00000000e+00, -0.00000000e+00, 0.00000000e+00, -0.00000000e+00, 5.18078151e-01, -1.61121546e-01, 0.00000000e+00, -4.24372747e-02, -5.16195159e-01, -4.26033726e-01, 2.12103848e-01, -5.75149604e-01, 3.74467667e-01, 1.51477181e-01, -0.00000000e+00, -0.00000000e+00, -2.59401414e-03, -6.81035927e-01, -1.35884843e-01, 1.20519287e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.19710887e-01, 0.00000000e+00, -0.00000000e+00, -0.00000000e+00, 4.32392321e-01, 5.95061672e-01, -0.00000000e+00, -7.34775902e-01, 3.92426845e-01, 0.00000000e+00, -0.00000000e+00, -1.29572606e-01, -3.02331023e-02, 1.14892681e-01])

このように、係数の値しか入っていません。そこで、前回sklearn.preprocessing.PolynomialFeaturesの仕様を調査したときに作成した関数を手直しして、この配列に対応する特徴量名を付けてみましょう

# 1番目からn番目までの素数をnp.arrayで返す関数
import sympy
def primes(n=0):
    ret_primes = list()
    for i in np.arange(n):
        ret_primes.append(sympy.prime(i+1))
    return ret_primes

# コラム名リストからコラムに対してユニークな素数が割り付けられたデータフレームを返す関数
def generate_df_prime_from_columns(columns_original = np.array(["a","b","c"])):
    data_original = np.array(primes(len(columns_original)))
    return pd.DataFrame(data=data_original[np.newaxis,:],columns=columns_original,index=["original"])

# PolynomialFeatures()の出力に対応したカラムを返す関数
def get_columns_for_PolynomialFeatures(poly=PolynomialFeatures(2),columns_from=["a","b","c"],power=False):
    df_from = generate_df_prime_from_columns(columns_from)
    columns_from = df_from.columns
    data_from = df_from.values
    data_poly = poly.fit_transform(df_from)
    # columnsをもう一度作り直す
    columns_poly = list()
    for i, data_poly_0 in enumerate(data_poly[0]):
        if (data_poly_0 == 1):
            columns_poly.append("bias")
        else:
            prime_dict=sympy.factorint(data_poly_0)
            keys = list(prime_dict.keys())
            columns_str = ""
            if power:
                # 累乗で書ける部分は累乗で書く(例:a^2)
                for j, key in enumerate(keys):
                    columns_str += columns_from[list(data_from[0]).index(key)]
                    if prime_dict[key] > 1:
                        columns_str += "^" + str(prime_dict[key])
                    if (j < len(keys)-1):
                        columns_str += "×"
            else:
                # 単純に×で項目をつなげていく(例:a×a×b)
                for j, key in enumerate(keys):
                    for k in np.arange(prime_dict[key]):
                        columns_str += columns_from[list(data_from[0]).index(key)]
                        if (j < len(keys)-1) | (k < prime_dict[key]-1):
                            columns_str += "×"
            columns_poly.append(columns_str)
    return columns_poly

前回はinvestigate_PolynomialFeatures という関数でしたが、今回は少し手直ししてget_columns_for_PolynomialFeaturesという名前にし、特徴量の名前(list形式)を返す関数にしています。

さて、get_columns_for_PolynomialFeaturesを使用して、係数の大きな順番にランキング形式で表示してみましょう。

# 係数を出力
best_estimator.named_steps["en"].coef_

columns_for_coef_ = get_columns_for_PolynomialFeatures(poly=PolynomialFeatures(best_estimator.named_steps["pf"].degree),
                                                       columns_from=df_features_.columns)
df_coef = pd.DataFrame(data=best_estimator.named_steps["en"].coef_,index=columns_for_coef_,columns=["coef"])
df_coef["abs coef"] = abs(df_coef["coef"])

# 係数の絶対値が0.1以上のものを絶対値の降順で表示
display(df_coef[df_coef["abs coef"]>0.1].sort_values('abs coef', ascending=False))

	coef	abs coef
bias	17.580475	17.580475
RM	2.504036	2.504036
LSTAT	-2.486765	2.486765
CRIM	-2.045977	2.045977
TAX	-0.950377	0.950377
INDUS	-0.923279	0.923279
RAD×TAX	0.919711	0.919711
AGE	-0.775039	0.775039
DIS	-0.767196	0.767196
TAX×LSTAT	-0.734776	0.734776
AGE×B	-0.681036	0.681036
TAX×PTRATIO	0.595062	0.595062
B	0.594885	0.594885
PTRATIO	-0.591576	0.591576
RM×LSTAT	-0.575150	0.575150
INDUS×INDUS	0.559900	0.559900
RM×RM	0.518078	0.518078
RM×TAX	-0.516195	0.516195
ZN	-0.486938	0.486938
CHAS×RM	-0.464449	0.464449
NOX	-0.454374	0.454374
ZN×PTRATIO	0.446274	0.446274
NOX×RM	-0.438023	0.438023
TAX×TAX	0.432392	0.432392
RM×PTRATIO	-0.426034	0.426034
NOX×DIS	0.422305	0.422305
CHAS×DIS	0.409789	0.409789
PTRATIO×PTRATIO	0.392427	0.392427
AGE×AGE	0.374468	0.374468
RM×B	0.212104	0.212104
CHAS×NOX	-0.204649	0.204649
CHAS×CHAS	0.170832	0.170832
RM×AGE	-0.161122	0.161122
ZN×ZN	0.154256	0.154256
AGE×DIS	0.151477	0.151477
AGE×LSTAT	-0.135885	0.135885
B×B	-0.129573	0.129573
DIS×DIS	0.120519	0.120519
CRIM×CRIM	0.118679	0.118679
LSTAT×LSTAT	0.114893	0.114893
ZN×RM	0.110348	0.110348
CRIM×B	-0.102504	0.102504

このように、どの特徴量の係数が大きいかが分かるようになりました。Elasticnetで回帰する前にStandardScalerで標準化しているので、基本的には係数が大きい特徴量が住宅価格を決めるのに重要となっていると考えられます。

さて、今回はBoston Housingデータセットを例に、sklearn.preprocessing.PolynomialFeaturesで生成した特徴量の係数名を得てみました。何かのお役に立てば幸いです。

*1:Boston Housingデータセットについては下記サイトを参考にしました。

「Boston Housing Prices - Evaluation & Validation」https://www.kaggle.com/samratp/boston-housing-prices-evaluation-validation
「Linear Regression on Boston Housing Dataset」https://towardsdatascience.com/linear-regression-on-boston-housing-dataset-f409b7e4a155
「機械学習とBoston Housing Data（3）」
https://kiidax.wordpress.com/2016/09/14/機械学習とboston-housing-data（3）/

*2:グリッドサーチはこちらのサイトを参考にいたしました。

「GridSearchCV(scikit-learn)によるチューニング」http://starpentagon.net/analytics/scikit_learn_grid_search_cv/

2019-08-01

sklearn.preprocessing.PolynomialFeatures の仕様を調査する

python 機械学習 scikitlearn

こんにちは、company-rapperです

初ポストです

sklearn.preprocessing.PolynomialFeaturesは簡単に特徴量同士の積や自分自身の累乗を作り出すことができて便利ですが、

そもそもどの順番で何の積を出しているんだろう!?

と思ったことはありませんか!?

本家本元のドキュメントを読んでも・・・

scikit-learn.org

いまいちどの順番で何を出すかの仕様がわかりません。

ネットで調べても、はっきりと説明してくれているところが私が調べた限りでは見つからず、

それならば自分がやってみようと思い立った次第です。

調査のカラクリは単純で、ユニークな素数をnumpy行列の各要素に割り当てて、PolynomialFeaturesで変換後のnumpy行列を得て、その項目ごとに素因数分解をして、「何と何をかけたものであるか」を調べます。

コード(Jupyter notebook形式)はGithubにアップしています。

github.com

結論としては、

バイアス→係数1個→係数2個→・・・の順
係数2個の中では、a×(a→b→c)→b×(b→c)→c×cの順
interaction_only = True (デフォルトはFalse)の時は同じ係数が2個以上登場しないもののみとなる

です。以降は、この結論に至るコードについて説明します!

まずは、必要なライブラリをインポートします

import pandas as pd
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
import sympy

ここでは、素数を出力したり、素因数分解ができるメソッドを備えたsympyを使用します。

調査のために、a,b,cというカラム名を持ったpandasのDataFrameを作成します。

# 1番目からn番目までの素数をnp.arrayで返す関数
def primes(n=0):
    ret_primes = list()
    for i in np.arange(n):
        ret_primes.append(sympy.prime(i+1))
    return ret_primes

# コラム名リストからカラムに対してユニークな素数が割り付けられたデータフレームを返す関数
def generate_df_prime_from_column(columns_original = np.array(["a","b","c"])):
    data_original = np.array(primes(len(columns_original)))
    return pd.DataFrame(data=data_original[np.newaxis,:],columns=columns_original,index=["original"])

# テスト
display(generate_df_prime_from_column())

テストの出力は

	a	b	c
original	2	3	5

となります。データとして1行3列で要素として、ユニークな素数を与えています。

sklearn.preprocessing.PolynomialFeaturesにこのDataFrameを入力すると、ユニークな素数同士の掛け算の結果がnumpy配列(ndarray)形式で返ってきます。このndarray形式の配列の要素を一つずつ素因数分解すると、何の項目同士を掛けたかが分かるという仕組みです:

def investigate_PolynomialFeatures(poly=PolynomialFeatures(2),columns_from=["a","b","c"],power=False):
    df_from = generate_df_prime_from_column(columns_from)
    columns_from = df_from.columns
    data_from = df_from.values
    data_poly = poly.fit_transform(df_from)
    # columnをもう一度作り直す
    columns_poly = list()
    for i in np.arange(data_poly.shape[1]):
        if (data_poly[0][i] == 1):
            columns_poly.append("bias")
        else:
            prime_dict=sympy.factorint(data_poly[0][i])
            keys = list(prime_dict.keys())
            column_str = ""
            if power:
                # 累乗で書ける部分は累乗で書く(例:a^2)
                for j in np.arange(len(keys)):
                    column_str += columns_from[list(data_from[0]).index(keys[j])]
                    if prime_dict[keys[j]] > 1:
                        column_str += "^" + str(prime_dict[keys[j]])
                    if (j < len(keys)-1):
                        column_str += "×"
            else:
                # 単純に×で項目をつなげていく(例:a×a×b)
                for j in np.arange(len(keys)):
                    for k in np.arange(prime_dict[keys[j]]):
                        column_str += columns_from[list(data_from[0]).index(keys[j])]
                        if (j < len(keys)-1) | (k < prime_dict[keys[j]]-1):
                            column_str += "×"
            columns_poly.append(column_str)
    return pd.DataFrame(data=data_poly,columns=columns_poly,index=["poly"])

sympy.factorint()は整数を与えると、素因数分解した結果をdict形式で返してくれます。例えば8なら、{2: 3} (2の3乗)と返してくれます。これを基に素数に対応する元のDataFrameのカラム名を取得して何の掛け算で構成されているかを再構成し、カラム名を作成しています。

さて、実際に動かしてみると、以下のような結果になります

degree=2のとき
	bias	a	b	c	a×a	a×b	a×c	b×b	b×c	c×c
poly	1.0	2.0	3.0	5.0	4.0	6.0	10.0	9.0	15.0	25.0

degree=2, interaction_only=Trueのとき
	bias	a	b	c	a×b	a×c	b×c
poly	1.0	2.0	3.0	5.0	6.0	10.0	15.0

degree=2, interaction_only=True, include_bias=Falseのとき
	a	b	c	a×b	a×c	b×c
poly	2.0	3.0	5.0	6.0	10.0	15.0

degree=3, interaction_only=False, include_bias=Falseのとき
	a	b	c	a^2	a×b	a×c	b^2	b×c	c^2	a^3	a^2×b	a^2×c	a×b^2	a×b×c	a×c^2	b^3	b^2×c	b×c^2	c^3
poly	2.0	3.0	5.0	4.0	6.0	10.0	9.0	15.0	25.0	8.0	12.0	20.0	18.0	30.0	50.0	27.0	45.0	75.0	125.0

degree=3, interaction_only=True, include_bias=Falseのとき
	a	b	c	a×b	a×c	b×c	a×b×c
poly	2.0	3.0	5.0	6.0	10.0	15.0	30.0

これをまとめると、

バイアス→係数1個→係数2個→・・・の順
係数2個の中では、a×(a→b→c)→b×(b→c)→c×cの順
interaction_only = True (デフォルトはFalse)の時は同じ係数が2個以上登場しないもののみとなる

ということが分かりました。ルールが明らかになれば、大したことじゃなかったんだなと思ってしまいます・・・

さて、このinvestigate_PolynomialFearures()ですが、単にsklearn.preprocessing.PolynomialFeaturesの挙動を調査するだけでなく、PolynomialFeatures適用後のcolumnが得られる優れものにもなるのです!

ということで、次回は調査目的ではなく、実用目的にフォーカスして少しコードを書き換えたもので、何か分析してみようと思います。

(次の投稿はいつになることやら・・・)

company-rapperでした!