統計 - 静かなる名辞

【python】matplotlibのboxplotで外れ値を表示しないようにする

Mon, 07 Oct 2019 20:24:16 +0900

はじめに

　matplotlibのboxplotを使うと簡単に箱ひげ図が描けます。ただし、デフォルト設定では外れ値が黒い円で表示されます。

　どんなデータでも、サンプル数が多いと一定数の外れ値は出てしまいます。ただ、図を見る人は気にするところですし、外れ値がたくさんあると見た目にも悪いので、何らかの処置が必要です。

　外れ値が描画されてしまうプログラムの例

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
x = np.random.normal(size=3*10**3)  # サンプル数に注目
plt.boxplot(x)
plt.savefig("result1.png")

result1.png　外れ値がたくさん出ている

　ということで、外れ値を表示させない方法を解説します。なお、この記事の内容は公式リファレンスに基づいています。

matplotlib.pyplot.boxplot — Matplotlib 3.1.1 documentation

シンプルに表示させない

　boxplotのキーワード引数のsymを使うと外れ値を「描画させない」設定が可能になります。具体的には、空文字列を指定します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
x = np.random.normal(size=3*10**3)
plt.boxplot(x, sym="")  # 変更点
plt.savefig("result2a.png")

result2a.png　外れ値が表示されなくなった

　注意点としては、これはあくまでも外れ値の描画を行わないだけで、外れ値そのものはデフォルト通り計算されるということが挙げられます。つまり、本来の最大値・最小値よりひげが短くなります。

　プロット上でだけ「消している」ということですね。

　参考：
4-3. 外れ値検出のある箱ひげ図 | 統計学の時間 | 統計WEB

　また、外れ値のマーカーに好きな記号を指定したりすることもできます。本来はこの用途で用いる引数です。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
x = np.random.normal(size=3*10**3)
plt.boxplot(x, sym="+")  # 変更点
plt.savefig("result2b.png")

result2b.png　マーカーを十字にしてみた

　使えるマーカーの一覧はこちらを参照。
matplotlib.markers — Matplotlib 3.1.1 documentation

外れ値の計算そのものをやめる

　（外れ値も含めた）本来の最大値・最小値に基づいてひげを出す場合は、whis="range"を指定します。

　そもそも外れ値の概念なしで箱ひげ図を描くというやり方になります。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
x = np.random.normal(size=3*10**3)
plt.boxplot(x, whis="range")  # 変更点
plt.savefig("result3.png")

result3.png　上の方法で非表示にするよりひげが長くなる

　この方が誤解を招く恐れがないので（つまり、外れ値検出をやっておいて外れ値をプロットしないという図は少しイレギュラーで伝わりづらい気がするので）、これでやるのがおすすめです。ただし、外れ値が少なくて、そんな極端に外れていないとき向きのやり方です（やたらひげが長いのも不格好ですから……）。

まとめ

　二つやり方がありますが、意味が微妙に違うし、実際にそれぞれで異なった結果になるので注意してください。

【python】相関係数行列をstatsmodelsを使って描く

Sun, 21 Jul 2019 01:54:42 +0900

はじめに

　相関係数行列を描く方法としては、pandasとseabornを使う方法などが一般的です。しかし、statsmodelsで行う方法も実は存在します。

pandas+seabornでやる場合

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(0)
df = pd.DataFrame({"A":np.random.randint(0, 10, size=(10,)),
                   "B":np.random.randint(0, 10, size=(10,)),
                   "C":np.random.randint(0, 10, size=(10,))})

df_corr = df.corr()
sns.heatmap(df_corr, vmax=1, vmin=-1, center=0)
plt.savefig("pandas_cm.png")

pandas_cm.png

　参考：pandas.DataFrameの各列間の相関係数を算出、ヒートマップで可視化 | note.nkmk.me

　というか相関係数の記事を書くときにstatsmodelsのリファレンスを検索していたらこれを見つけたので、試してみようと思った次第です。どれくらい使い物になるのでしょうか。

使い方

　ずばり、これです。

statsmodels.graphics.correlation.plot_corr — statsmodels v0.10.1 documentation

　名前が覚えづらいのが最大の難点で、他は普通に使えます。というか、seaborn.heatmapとそんなに変わりません。

　公式で紹介されている使い方。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import statsmodels.graphics.api as smg
>>> hie_data = sm.datasets.randhie.load_pandas()
>>> corr_matrix = np.corrcoef(hie_data.data.T)
>>> smg.plot_corr(corr_matrix, xnames=hie_data.names)
>>> plt.show()

　相関行列の計算からやってくれるわけではなく、例ではnumpyで計算しているようです。

　これを踏まえて書いたコード。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.correlation import plot_corr

np.random.seed(0)
df = pd.DataFrame({"A":np.random.randint(0, 10, size=(10,)),
                   "B":np.random.randint(0, 10, size=(10,)),
                   "C":np.random.randint(0, 10, size=(10,))})

df_corr = df.corr()
plot_corr(df_corr, xnames="ABC")
plt.savefig("statsmodels_cm.png")

　満足の行く図にするためには、最低限xnamesを指定する必要があります。seaborn.heatmapでいうxticklabels, yticklabelsですね。iterableを渡せば良いようです。

　結果。

statsmodels_cm.png

　見た目が少し違いますね。デフォルトのカラーマップはこっちの方が良いかも。

　それだけ？　と思ってリファレンスを見直したら、面白そうな引数が2つだけありました。

title
str, optional
The figure title. If None, the default (‘Correlation Matrix’) is used. If title='', then no title is added.

normcolor
bool or tuple of scalars, optional
If False (default), then the color coding range corresponds to the range of dcorr. If True, then the color range is normalized to (-1, 1). If this is a tuple of two numbers, then they define the range for the color bar.

　normcolorはカラーマップの範囲を-1～1に一致させてくれるので、その方が良いと思います。同じことをseabornでできないか調べましたが、ぱっと見なさそうなので、アドバンテージです。

　両方設定してみましょう。
　

plot_corr(df_corr, xnames="ABC", title="random", normcolor=True)

statsmodels_cm.png

　それなりに満足感のある結果。

まとめ

　どうということはありませんが、欠点もなさそうだし、使いやすいのでこちらでやってもいいはず。

　難点は、やはり名前が覚えづらいこと（importを書くときに迷う）くらいでしょうか。

pythonで相関係数を計算する方法いろいろ3種類

Thu, 18 Jul 2019 04:16:46 +0900

はじめに

　pythonで相関係数を計算する方法はいろいろあります。確認したら、主要ライブラリだけで3つありました。

　いろいろあるということは用途によって使い分けられるということなので、淡々と書いていきます。

　なお、念のために断っておくと、ここで書いている「相関係数」はすべて「ピアソンの積立相関係数」です。順位相関などはまた別に調べてください（ただしpandasを使う方法だと出せます）。

はじめに
データの確認
numpyでやる
pandasでやる
scipyを使う
あと思ったこととか
まとめ

データの確認

　予め以下のようなデータを定義しておきます。

>>> import numpy as np
>>> np.random.seed(0)
>>> x = np.arange(0, 10, 0.1)
>>> y = x + np.random.normal(size=x.shape)

　散布図にプロットして確認。

>>> import matplotlib.pyplot as plt
>>> plt.scatter(x, y)
<matplotlib.collections.PathCollection object at 0x7f31aa415f28>
>>> plt.savefig("fig.png")

fig.png

　もう少しサンプル数が少なくても良かったような気もしますが、せっかく定義したのでこれでやります。

numpyでやる

　numpyの場合はnp.corrcoefで相関係数「行列」を出してくれます。

>>> np.corrcoef(x, y)
array([[1.        , 0.94129622],
       [0.94129622, 1.        ]])

　0.9以上なので強い相関があるみたいです。「行列」が出てくるので、単に相関係数がほしいときは適当に取り出します。

>>> np.corrcoef(x, y)[0, 1]
0.9412962237004372

　あまりスマートではないので、本当に相関係数「行列」がほしいときに使います。

numpy.corrcoef — NumPy v1.17 Manual

pandasでやる

　pandasでもnumpyと同じことができるようです。

>>> import pandas as pd
>>> df = pd.DataFrame({"x":x, "y":y})
>>> df.corr()
          x         y
x  1.000000  0.941296
y  0.941296  1.000000

　行と列に名前がついて使いやすくなったと思います。また、ピアソン以外の相関係数も、kendall, spearmanをmethod引数に渡すことができ、なんならcallableで任意の関数で計算することもできるといった使いやすさがあります。多機能ですね。

pandas.DataFrame.corr — pandas 0.25.1 documentation

　あと、相関係数「行列」がほしいときはpandasを経由した方が便利でしょうか。seabornに投げて可視化するときに、行・列の名前を考慮してくれるので、便利そうです。

pandas.DataFrameの各列間の相関係数を算出、ヒートマップで可視化 | note.nkmk.me

scipyを使う

　漢は黙ってscipy、という価値観が私にはあります。

>>> from scipy import stats
>>> stats.pearsonr(x, y)
(0.941296223700437, 5.153124094421605e-48)

　勝手に両側検定をやってp値を出してくれています（結果のtupleの0から数えて1つめ）。

scipy.stats.pearsonr — SciPy v1.3.0 Reference Guide

　検定やってくれるのはいいですね。普通は別途やる必要があると思います。

あと思ったこととか

なんで標準のstatisticsで用意されてないの
statsmodelsは高度な機能はいろいろ提供しているくせに、ただの相関係数の出し方がいくらググっても出てこないのはなんで。リファレンスすごく読みづらいし。あるかもしれないけど諦めた

まとめ

　まあ、3つあればいいか……行列がほしいときは楽そうなのはpandas、単に数字がほしければscipyという使い分けになりそうですね。

scikit-learnで目的変数を対数変換したりするTransformedTargetRegressor

Sun, 14 Jul 2019 04:33:24 +0900

はじめに

　経済系の分析などで、目的変数を対数変換して分析するというケースがあります。scikit-learnはそのようなケースもサポートしています。

　どうやったらいいのかわからなくて、自分で変数を変換している人も中にはいるかと思いますが、モデル構築まではなんとかなっても、予測のことまで考えると不便になります。うまくPipelineなどで自動化できるといいのですが、普通のやり方では目的変数は処理してくれません。しかし、TransformedTargetRegressorなら大丈夫です。

目的変数の対数変換

　sklearn.compose.TransformedTargetRegressorを使います。ほとんど紹介を見かけないのですが、実際これでできます。

>>> import numpy as np
>>> from sklearn.linear_model import LinearRegression
>>> from sklearn.compose import TransformedTargetRegressor
>>> tt = TransformedTargetRegressor(regressor=LinearRegression(),
...                                 func=np.log, inverse_func=np.exp)
>>> X = np.arange(4).reshape(-1, 1)
>>> y = np.exp(2 * X).ravel()
>>> tt.fit(X, y) 
TransformedTargetRegressor(...)
>>> tt.score(X, y)
1.0
>>> tt.regressor_.coef_
array([2.])

sklearn.compose.TransformedTargetRegressor — scikit-learn 0.21.2 documentation

　なかなか特殊な感じですね。モデルを作るときにfuncとinverse_funcを渡すのがミソで、学習はfuncで変換された目的変数に対して行われます。予測のときはinverse_funcで逆変換されて出てくるので、余計な手間をすべて省くことができます。

やってみる

　ということで、対数変換した方がうまくいくようなデータを作って線形回帰してみます。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.compose import TransformedTargetRegressor

def main():
    np.random.seed(0)
    x = np.linspace(0, 5, 100)
    y = 20 + np.exp(x + 
                    np.random.normal(scale=0.5, size=x.shape))
    X = x.reshape(-1, 1)

    lr = LinearRegression()
    lr_logy = TransformedTargetRegressor(
        regressor=lr, func=np.log, inverse_func=np.exp)

    plt.scatter(x, y, c="b", alpha=0.1)
    for name, model in [("Linear", lr), 
                        ("Linear-logy", lr_logy)]:
        model.fit(X, y)
        pred = model.predict(X)
        plt.plot(x, pred, label=name)
    plt.legend()
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

result.png

　非等分散だったりして条件が悪いのですが、線形回帰で素直にやるのと比べると良いフィッティングを見せています。

おまけ：説明変数を対数変換したいとき

　FunctionTransformerを使ってください。

sklearn.preprocessing.FunctionTransformer — scikit-learn 0.21.2 documentation

>>> import numpy as np
>>> from sklearn.preprocessing import FunctionTransformer
>>> transformer = FunctionTransformer(np.log1p, validate=True)
>>> X = np.array([[0, 1], [2, 3]])
>>> transformer.transform(X)
array([[0.        , 0.69314718],
       [1.09861229, 1.38629436]])

5.3. Preprocessing data — scikit-learn 0.21.2 documentation

　Pipelineで組み合わせれば、説明変数・目的変数ともに対数変換でloglogというのも簡単です。

まとめ

　対数変換した方がうまくいくようなデータのときは試してみましょう。もちろん対数以外の変換でもいけるので、目的変数を変換したいときはこれでいいと思います。

ロジスティック回帰が線形分離不可能な分類問題を解けないことの説明

Sun, 07 Jul 2019 03:36:21 +0900

はじめに

　ロジスティック回帰が線形分離不可能な分類問題を解けないことは有名な話です。だけど、「いや解けるだろ」「なんで解けないの？？？」と言われてしまうことがあるので*1、それができないことを説明しておこうと思います。

　なお、この記事はこちらの記事を参考にしています。

誤った図解から学ぶロジスティック回帰の性質 - ill-identified diary

　書きたいことは言い尽くされている感もあるので、こういう結論に至る過程を数式で書きます*2。

y=0.5を代入すればいい

　さて、説明変数、目的変数、パラメータ、などを適当に定めたとします。すると、ロジスティック回帰の予測式はこんなやつになります（は適当に学習できたとする）。

\begin{align}
\hat{y} = \frac{1}{1 +
\mathrm{e}^{-(\beta_0 + \sum_{i=1}^{n}\beta_i x_i)}
}
\end{align}

　書き方の流儀はいろいろあると思いますが（を使うとか）、今回は上の式で行きます。

　さて、今回は分離超平面の式に興味があるのでしたね。分離超平面ってどこ？　というと、のところです。なので、淡々と式を書き換えます。

\begin{align}
0.5 = \frac{1}{1 +
\mathrm{e}^{-(\beta_0 + \sum_{i=1}^{n}\beta_i x_i)}
}
\end{align}
　
　両辺を逆数にします。

\begin{align}
2 = 1 + \mathrm{e}^{-(\beta_0 + \sum_{i=1}^{n}\beta_i x_i)}
\end{align}

　とりあえず邪魔な1を反対側に移す。

\begin{align}
1 = \mathrm{e}^{-(\beta_0 + \sum_{i=1}^{n}\beta_i x_i)}
\end{align}

　両辺の対数を取る。

\begin{align}
0 = -(\beta_0 + \sum_{i=1}^{n}\beta_i x_i)
\end{align}

　マイナスは必要ないので消しましょう。

\begin{align}
0 = \beta_0 + \sum_{i=1}^{n}\beta_i x_i
\end{align}

　もうだいたい終わってる気もしますが、たとえばとして適当に式を変形します。

\begin{align}
\beta_0 + \beta_1 x_1 + \beta_2 x_2 = 0
\end{align}

　……はい、これは「直線の式」ですね。

　が増えると係数と変数が増えていきますが、いずれにせよ線形の式なのは間違いありません。とかみたいなのが出てくる余地はありません。

　「分離超平面」がかのような式で表わせる以上、線形分離不可能な分類問題は解けない、ということです。

非線形の問題も解く方法

　あくまでも「分離境界が線形にならないといけない」というだけなので、データを非線形変換して解けるような空間に写像すればできます。代表的な方法は多項式を使うことです（SVMの多項式カーネルなんかと同じですが、明示的に特徴量空間を計算するのが相違点です）。

　ということで、こちらの記事を御覧ください。どれくらい非線形でも行けるのかが書いてあります。

非線形がなんだ！ロジスティック回帰＋多項式でやってやる！ - 静かなる名辞

*1:……

*2:誰でも納得するから

コサイン距離は距離じゃないんだから、勘違いしないでよねっ！

Fri, 05 Jul 2019 03:09:35 +0900

き、記事タイトルに意味なんてないんだからねっ！

　自然言語処理などでお馴染みのコサイン類似度。これを1から引いたものを「コサイン距離」と称している文献も散見されますが、この「コサイン距離」は距離としての性質を満たしません。

　それがどういうことなのかをこの記事で説明していきます。

コサイン類似度のことくらい自分で調べなさいっ！

　まず前提となるコサイン類似度については、親切に解説しているサイトが他にたくさんあるので、そちらに譲ります。

　たとえばここなどがいいでしょう。

コサイン類似度

　コサイン類似度はベクトル同士の類似度であり、要するに単なる内積（をノルムで正規化したもの）です。これは-1から1の区間を取ります。1なら「最も似ている（同じベクトル）」、-1なら「最も似ていない（反対向き）」という性質を持ちます。

　これを1から引くことで、0なら「最も似ている」、2なら「最も似ていない」に変換したものが「コサイン距離」です。

距離の定義を知らないの？　しょ、しょうがないから教えてあげるわ

　さて、距離という言葉というか概念は実は数学的にちゃんと定義できます。かいつまんで書くと、関数が以下の条件（距離の公理といいます）を満たすとき、その関数を距離関数あるいは距離と言えます。

\begin{align}
d(x,y) &>& 0\\
x&=&y\Leftrightarrow d(x, y) = 0\\
d(x, y) &=& d(y, x)\\
d(x, z) &\leq& d(x, y) + d(y, z)
\end{align}

参考
距離空間 - Wikipedia
第６回距離の公理：ねむねこ幻想郷：So-netブログ
 距離とは (キョリとは) [単語記事] - ニコニコ大百科

　数式で見ると難しく見えるかもしれませんが、この式はそれぞれ

距離は負にはならない（非負性）
同じ点同士の距離は0、距離が0の点は同じ点
x,yの間の距離について、距離を測る起点を逆にしても距離は変わらない（対称性）
x,zとまっすぐ行くときと比べて、yに寄り道すると必ずトータルの経路は長くなる（三角不等式）

　ということを言っているだけなので、概念的には簡単です。

　こういうものを満たすと距離と呼べる、ということですね。

　「コサイン距離」はどれを満たさないのでしょうか？

わからないの？　……ばか

　「コサイン距離」は2番目のと、4番目の三角不等式を満たしません。

　2番目を説明するのは簡単で、元のコサイン類似度はベクトル間の角度にしか興味を持たない性質があります。なので、たとえば二次元ベクトルととか、とは同じ距離になります。

　4番目については、反例を挙げてみましょう。

すごく単純な例

　特に凝ったことはしていません。この図において、単純なユークリッド距離を考えると、

A-B間, B-C間の距離：1
A-C間の距離：

　となり、こういうのが三角不等式を満たしている場合です。A-B-Cとたどる経路の長さは2になるので、A-Cとたどるより長い距離をたどることになります。

　では、「コサイン距離」では？　というと、

A-B間, B-C間の「コサイン距離」：約0.293
A-C間の距離：1

　となり、A-B-Cとたどることで約0.586になりますからA-Cと直接たどるより短い距離で行けてしまうことになります。つまり、三角不等式を満たさないので、「コサイン距離」は距離ではないということになります。

距離として扱うと困るのかって？　……困るに決まってるじゃないっ、わからずや！

　データ分析などで、距離を使うことを前提としている手法で「コサイン距離」を使うと、不都合なことが起きる可能性があります。

　みんなが大好きなirisのデータを多次元尺度構成法、MDSで可視化してみましょう。Pythonで書くとこうなります。

import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform
from sklearn.datasets import load_iris
from sklearn.manifold import MDS

def main():
    iris = load_iris()

    A = squareform(pdist(iris.data, "euclidean"))
    mds = MDS(n_components=2, dissimilarity="precomputed",
              n_init=10, max_iter=500)
    X_2d = mds.fit_transform(A)

    for i, target in enumerate(iris.target_names):
        mask = iris.target == i
        plt.scatter(X_2d[mask,0], X_2d[mask,1], label=target)

    plt.xlim(X_2d[:,0].min() - 1, X_2d[:,0].max() + 1)
    plt.ylim(X_2d[:,1].min() - 1, X_2d[:,1].max() + 1)
    plt.title("MDS stress:{:.4f}".format(mds.stress_))
    plt.legend()
    plt.savefig("iris_euclidean_mds.png")
   
if __name__ == "__main__":
    main()

iris_euclidean_mds.png

　そんなに申し分はなさそうな結果ですね。

　「コサイン距離」でもやってみます。といっても、親切なことにscipyが「コサイン距離」を標準でサポートしているので、

    A = squareform(pdist(iris.data, "cosine"))

　とすれば一発でできます。あとはせいぜい出力ファイル名を変えておきます。
（plt.savefig("iris_cosine_mds.png")としました。）

scipy.spatial.distance.pdist — SciPy v1.3.0 Reference Guide

iris_cosine_mds.png

　なんかよくわからないことになりました。念のために中心付近にズームしてみます（plt.xlimとplt.ylimで調整）。

iris_cosine_mds_zoomed.png

　考えてみれば当然の結果で、コサイン類似度は-1から1のレンジを取ります。ということは、「コサイン距離」の最大値は2にしかならないのです。なので、遠い点が表現できなくなり、とても小さい範囲に押し込められます。

　また、「コサイン距離」では向きが同じで長さの違うベクトル同士を区別できません。昔作ったirisの主成分分析のバイプロットを持ってくると、

irisのバイプロット

【python】pythonで主成分分析のバイプロット - 静かなる名辞

　グループ間の差異は概ね第一主成分に、グループ内での差異は第二主成分にあらわれています。そして、第一主成分とほぼ同じ方向を向いている2つの変数、そうでもない2つの変数があることがわかります。

　品種が違うと各変数の相対的な比率が変わる反面、同じ品種同士では各変数の相対的な比率はさほど変わらない（全体的に大きかったり小さかったりという個体差があるだけ）と想定すれば、結果が一直線上に並ぶのもなんとなく理解できる気がします。

「じゃあどう呼べば良いのか」って？　そんなの自分で考えてよね！

　「コサイン距離」に変わる呼称方法ですが……

　ま、常識的に考えると、コサイン非類似度でいいのではないでしょうか。

わかったなら感謝しなさい。……え、ありがとう？　べ、べつに喜ばれても嬉しくなかんないんだからっ！

　安易に「コサイン距離」という言葉を使ってはいけないこと、また、距離として扱うと問題になるというか、イマイチな結果を招く可能性があることがこの記事でわかっていただけたら、嬉しいです。

　あと、ツンデレ風の章タイトルにしたことに対して今更ながら後悔の念を感じ始めているのですが（自分で見返してもかなり痛い）、下書きに放り込んで一晩寝たらたぶん投稿する勇気がなくなっていると思うので、蛮勇を奮ってこのまま後悔公開することにします。

【python】ロジスティック回帰で確率値で学習させる

Sun, 30 Jun 2019 22:04:06 +0900

はじめに

　ロジスティック回帰は回帰という名前なのにほとんど二項判別に使われますが、たまに本当に回帰に使うときもあります。0.1とか0.4とか0.6のような目的変数を使ってモデルを作る、というケースです。

　ちょっとした目的で必要になるかもしれないと思ってやろうとしたら、意外と手間取ったのでメモしておきます。

データ

　たとえば「普及率」のようなデータに対してあてはめを行うとき、こういうケースが出てきます。

　こちらで紹介されている、日本のカラーテレビ普及率のデータを使います。

データ解析・マイニングとR言語

　説明変数が年、目的変数が普及率です。

　とりあえずこんな配列にしておきます。

import numpy as np

x = np.array([1966, 1967, 1968, 1969, 1970,
              1971, 1972, 1973, 1974, 1975,
              1976, 1977, 1978, 1979, 1980,
              1981, 1982, 1983, 1984]).reshape(-1, 1)

y = np.array([0.003, 0.016, 0.054, 0.139, 0.263,
              0.423, 0.611, 0.758, 0.859, 0.903,
              0.937, 0.954, 0.978, 0.978, 0.982,
              0.985, 0.989, 0.988, 0.992])

scikit-learnでは（たぶん）できない

　誰でもまっさきに思いつく方法は、sklearnのLogisticRegressionを使うことです。しかし、これは

Logistic Regression (aka logit, MaxEnt) classifier.
sklearn.linear_model.LogisticRegression — scikit-learn 0.21.2 documentation

　と書いてあるとおり、判別用のモデルです。ユーザガイドもひたすら判別の話をしているだけです。

1.1. Generalized Linear Models — scikit-learn 0.21.2 documentation

　まあでも、もしかしたらできるかもしれないので、やってみましょう。

import numpy as np
from sklearn.linear_model import LogisticRegression

x = np.array([1966, 1967, 1968, 1969, 1970,
              1971, 1972, 1973, 1974, 1975,
              1976, 1977, 1978, 1979, 1980,
              1981, 1982, 1983, 1984]).reshape(-1, 1)

y = np.array([0.003, 0.016, 0.054, 0.139, 0.263,
              0.423, 0.611, 0.758, 0.859, 0.903,
              0.937, 0.954, 0.978, 0.978, 0.982,
              0.985, 0.989, 0.988, 0.992])

lr = LogisticRegression()
lr.fit(x, y)  # => ValueError: Unknown label type: 'continuous'
lr.predict(x)

　しってた。

statsmodelsでやる

　仕方がないので、statsmodelsを使います。statsmodelsはPythonでR風のことをするためのライブラリなので、参考ページと同じことができるはずです。

　APIをぜんぜん把握していないので、qiitaの解説記事を見ながらやります。

Statsmodels でロジスティック回帰を行う際の注意点 - Qiita

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

x = np.array([1966, 1967, 1968, 1969, 1970,
              1971, 1972, 1973, 1974, 1975,
              1976, 1977, 1978, 1979, 1980,
              1981, 1982, 1983, 1984]).reshape(-1, 1)

y = np.array([0.003, 0.016, 0.054, 0.139, 0.263,
              0.423, 0.611, 0.758, 0.859, 0.903,
              0.937, 0.954, 0.978, 0.978, 0.982,
              0.985, 0.989, 0.988, 0.992])

x_c = sm.add_constant(x)  
# ↑interceptのためにやらないといけないらしい（えぇ…）
lr = sm.Logit(y, x_c)
lr_result = lr.fit()
print(lr_result.params)
print(lr_result.summary())

y_pred = lr.predict(lr_result.params)
plt.scatter(x.ravel(), y_pred, c="b", alpha=0.2)
plt.plot(x.ravel(), y_pred, c="b")
plt.savefig("result.png")

　y, xって書くのがきもいとか、どうでもいいところが気になります。

　これはこれでRに慣れてる人にはいいと思うのですが、scikit-learnライクなAPIも用意してくれていたらと思わなくはありません。

　結果

Optimization terminated successfully.
         Current function value: 0.180377
         Iterations 9
[-1.23730786e+03  6.27547565e-01]
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                   19
Model:                          Logit   Df Residuals:                       17
Method:                           MLE   Df Model:                            1
Date:                Sun, 30 Jun 2019   Pseudo R-squ.:                  0.7003
Time:                        22:00:31   Log-Likelihood:                -3.4272
converged:                       True   LL-Null:                       -11.435
Covariance Type:            nonrobust   LLR p-value:                 6.284e-05
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const      -1237.3079    591.413     -2.092      0.036   -2396.456     -78.159
x1             0.6275      0.300      2.092      0.036       0.040       1.215
==============================================================================

結果

　なんかよくわからないけど、一応できているみたいです。

まとめ

　いろいろ分析していくとstatsmodelsもけっきょく必要になるときがあるので、慣れた方が良いのかなぁと思ったりしました。

【python】statsmodelsでt検定する方法

Fri, 24 May 2019 00:58:37 +0900

はじめに

　この前はscipyでやる方法をまとめたわけですが、

【python】scipyでt検定する方法まとめ - 静かなる名辞

　片側のオプションほしいなーと思ったのでstatsmodelsに浮気することにしました。

使い方の概要

　対応のないt検定はこれです。

statsmodels.stats.weightstats.ttest_ind — statsmodels v0.10.1 documentation

　引数は以下のようなものです。

statsmodels.stats.weightstats.ttest_ind(x1, x2,
alternative='two-sided', usevar='pooled', weights=(None, None), value=0)

　x1とx2が検定対象です。オプションで重要そうなのはalternativeとusevarです。

alternative

　"two-sided"と"larger"と"smaller"が渡せます。後ろの2つは片側のオプションです。

　こうしちゃうと、どっちがどっちに対して大きい/小さいのかわからん、とか考えないのかしら（x1 larger/smaller than x2（x1はx2より大きい/小さい）です。）。

usevar

　"pooled"と"unequal"が渡せます。スチューデントのt検定とウェルチのt検定に対応しているはずです。

　難しくはないけど、いちいちstrで重要なオプションを渡さないといけないのは少し面倒ですかね。間違った値にするとエラー出してくれるみたいなので、そういう意味では親切ですが。

実際にやってみる

　完全に前回と同じノリでやります。

　分布をとりあえず作る。

>>> from scipy import stats
>>> from statsmodels.stats.weightstats import ttest_ind
>>> d1 = stats.norm(loc=0, scale=1)
>>> d2 = stats.norm(loc=1, scale=1)

　やる。

>>> a = d1.rvs(3)
>>> b = d2.rvs(3)
>>> a
array([ 0.90831032, -0.88621515,  0.09060862])
>>> b
array([1.87384828, 2.92258811, 0.8539749 ])
>>> ttest_ind(a, b)
(-2.333627132285731, 0.07993392828229898, 4.0)

　t統計量、p値、自由度が返るみたいですね。

　ここで仮に有意水準0.05とすると、もう少しでいけそうだったけど切れなかったという残念な例です（そうなるまで何回か回しました）。そこで、片側検定にしてみます。

>>> ttest_ind(a, b, alternative="smaller")
(-2.333627132285731, 0.03996696414114949, 4.0)

　教科書通り半分のp値になって、めでたく「有意差」が出ました。教科書には「こういうことはやるな」と書いてあると思います。気をつけましょう。

　なお、片側の方向を間違えて指定した場合は、

>>> ttest_ind(a, b, alternative="larger")
(-2.333627132285731, 0.9600330358588506, 4.0)

　だいたい大きいp値になるので気づきます。でも、元のp値が0.42とかだったりすると案外気づかないかもしれない（どっちにしろ有意差ではないし、実害ないかもですが・・・）。

まとめ

　こっちの方が高機能だし、これでいいんじゃ？　という感じもします。でもscipy入れててもstatsmodels入れてない人は多いと思うので、微妙っちゃ微妙ですね。

【python】scipyでt検定する方法まとめ

Fri, 17 May 2019 22:38:49 +0900

概要

　いっっっつも使い方を忘れて調べているので、自分で備忘録を書いておくことにしました。

概要
t検定の概要
使う関数
引数と注意点
やってみる
まとめ
statsmodelsでやる方法

t検定の概要

　t検定と一口に言っても実際はいろいろありますが、今回やるのは2群の標本の平均に差異があるかどうかの検定です。帰無仮説は「両者の平均に差はない」、対立仮説は「両者の平均に差がある」です。

　詳しいことはwikipediaとかを見てください（手抜き）。

t検定 - Wikipedia

使う関数

　scipyのt検定を行う関数としては、

scipy.stats.ttest_ind
scipy.stats.ttest_rel

　の2つがあります。ttest_indは対応のないt検定、ttest_relは対応のあるt検定で使えます。

　使い所が多いのは対応のないt検定を行うttest_indの方なので、こちらだけ取り扱います。

引数と注意点

　いろいろあります。

scipy.stats.ttest_ind(a, b, axis=0, equal_var=True, nan_policy='propagate')

　a,bは普通にデータの入った1次元配列を渡して使うことが多いでしょう。axisという引数があることから想像が付く通り、多次元配列でも渡せるようです。使ったことはありません。

　equal_var=Trueだとスチューデントのt検定、equal=var=Falseだとウェルチのt検定です。これは等分散かどうかに関わらずウェルチのt検定で良いという話題があるので、Falseを指定してやると良いと思います。

等分散か否かに関わらずウェルチの t 検定を使う (べきである) - Qiita

　他の引数はあまり重要ではないので、説明を省略します。

　結果は(t統計量, p値)というtupleっぽいオブジェクト*1で返ります。p値が設定した有意水準（たとえば0.05）より小さいときに有意差があったと言えます（不慣れだと毎回「どっちだっけ？」と思うポイント）。

　また気になる点として、t検定は母集団が正規分布に従うというけっこうきつい仮定を置いています。しかし、実はあまり気にする必要はないという議論もあります。

実際にはが正規分布でなくても，が大きければ中心極限定理によりは正規分布に近づくので，この検定は母集団が正規分布かどうかには鈍感です。データの分布が正規分布かどうかの検定をしてから t検定を行う必要はまったくありません。
t検定

　そういうことらしいです。

　なお、scipyのt検定に片側検定のオプションはありません。両側検定の結果から計算するか、他のライブラリ（statmodelsなど）でやることになります。両側検定の結果から計算する場合は、

t, p = stats.ttest_ind(male, female, equal_var=True)
pval3 = p
pval2 = p / 2.0
pval1 = 1.0 - pval2
if t < 0.0:
    w = pval2
    pval2 = pval1
    pval1 = w
pythonのpandasによる簡単な統計処理：第３回 F検定，t検定その他

　みたいな感じになるようです。

やってみる

　まず適当な確率分布のオブジェクトを生成する。

>>> from scipy import stats
>>> d1 = stats.norm(loc=0, scale=1)
>>> d2 = stats.norm(loc=1, scale=1)

　N(0, 1)とN(1, 1)です。

　標本はこんな感じで取れます。

>>> d1.rvs(10)
array([ 0.18711344,  0.3534579 , -0.52046706,  0.47855615, -0.51033227,
        0.70266909,  0.19253524,  0.28232341,  1.24373963, -0.70771188])

　参考：
scipyで確率分布のサンプルと確率密度関数を生成する - 静かなる名辞

　念のために正しいパラメータになっているか確かめます。

>>> np.mean(d1.rvs(1000))
0.031232377764520782
>>> np.var(d1.rvs(1000))
0.9921694887020086
>>> np.mean(d2.rvs(1000))
0.97389464006143
>>> np.var(d2.rvs(1000))
1.0268883977332324

　大丈夫そうなのでt検定します。有意水準0.05とします（先に決めるのがルールなので・・・）。

　最初は標本サイズ3でやってみます。

>>> a = d1.rvs(3)
>>> b = d2.rvs(3)
>>> a
array([-1.29621283,  0.42129238, -0.13701242])
>>> b
array([ 0.81419163,  1.21399486, -1.40737252])
>>> stats.ttest_ind(a, b, equal_var=False)
Ttest_indResult(statistic=-0.5672127490081906, pvalue=0.6064712381602329)

　pvalue=0.6064712381602329で、0.05より圧倒的におおきいので有意差なしということになります。サンプルサイズが少なすぎて有意差が出せないのです。

　10まで増やしてみます。

>>> c = d1.rvs(10)
>>> d = d2.rvs(10)
>>> stats.ttest_ind(c, d, equal_var=False)
Ttest_indResult(statistic=-2.8617115251996275, pvalue=0.011903486818782736)

　今度は出ました。ただし何回かやると有意になったりならなかったりするので、出方のばらつき次第で変わる可能性があります。

　サンプルサイズの見積もりは以下の方法があるそうです。

幾つデータが必要か？―平均値の差の検定 | ブログ | 統計WEB
　※引用者注：
　は標準偏差、は期待される二群間の差
　上式は有意水準5%の設定の場合に80%の検出力になるサンプル数
　（詳細はリンク先で読んでください）

　今回は標準偏差、二群間の差ともに1という簡単な設定なので、16サンプルあれば良い計算です。

>>> sum(stats.ttest_ind(
...       d1.rvs(16), d2.rvs(16), equal_var=False)[1] < 0.05
...       for _ in range(100))
78

　よさそうです。

まとめ

　簡単なことなのですが、割とやり方を忘れやすいので書き記しました。これで今後は忘れないで済むでしょう（フラグ）。

statsmodelsでやる方法

　statsmodelsでもできるので、こちらの記事も参考にしてください。片側検定が簡単にできるなどのメリットがあります。

【python】statsmodelsでt検定する方法 - 静かなる名辞

*1:厳密にはTtest_indResultという型だが、tupleの派生クラスで実質的にtupleとして扱えるので気にしなくて良い。中身はともにfloat

【python】scipyで階層型クラスタリングするときの知見まとめ

Mon, 11 Feb 2019 03:59:43 +0900

はじめに

　scipyの階層型クラスタリングを使う機会がありましたが、使い方がわかりづらいと思ったのでまとめておきます。

はじめに
関数がいっぱいある
使い方
- linkage
- fcluster
- cophenet
- dendrogram
実践編
- データを作る
- 手法を選ぶ
クラスタに分ける
- デンドログラムを描く
- 遊ぶ
まとめ

関数がいっぱいある

　いっぱいあるんですよ。

Hierarchical clustering (scipy.cluster.hierarchy) — SciPy v1.3.0 Reference Guide

　私の数え間違えがなければ31個。多いですね。

　とはいえ、本質的なもの（実際によく使うもの）は以下くらいです。

linkage

　実際に階層型クラスタリングを行う。これがないと始まらない。

fcluster

　任意の深さで木を切り、クラスタに分割する。

cophenet

　yを渡すとコーフェン相関係数なる評価指標を出してくれるらしい。

dendrogram

　デンドログラムを描画する。

　他に、こんな関数があります。

fclusterdata, leaders

　fclusterの処理と絡む便利関数。

single, complet , average,...など

　クラスタリングアルゴリズムに対応。すべてlinkageの引数で文字列を使って指定できるので実際にこれらの関数を使うことはない。

　他にも各種の数値を出せたり、MATLABのフォーマットと相互変換してくれたり、ポインタのリンクで表現された構造化データに変換するクラスだったり、いろいろなものが押し込まれています。必要に応じてリファレンスから探してくれば良いので、それぞれ述べることはしません。

使い方

　とりあえず上で挙げた実際によく使う関数を中心に説明していきますよ。

linkage

　データをlinkageに通すことで階層型クラスタリングが行えます。返り値として木の情報を表す配列が返ります。それに対して用意されている関数であれこれ処理していくというのが基本的な流れです。

scipy.cluster.hierarchy.linkage(y, method='single', metric='euclidean', optimal_ordering=False
scipy.cluster.hierarchy.linkage — SciPy v1.3.0 Reference Guide

　yはデータですが、基本的にはscipy.spatial.distance.pdistで作れる距離行列のフォーマット（一般的な正方距離行列とは異なるので注意）で渡してくれ、ということになっています。

scipy.spatial.distance.pdist — SciPy v1.3.0 Reference Guide

　でもshape=(サンプル数, ベクトル次元数)のnumpy配列も受け取ってくれるので、それほど気を配る必要はありません。

　気を配らないといけないのはむしろ距離行列を入れてクラスタリングしたい場合で、その場合は正方行列として入れようとすると正しく処理されません。pdistのフォーマットに変換する必要があります。これはscipy.spatial.distance.squareformで変換できるので、そうしてください。

scipy.spatial.distance.squareform — SciPy v1.3.0 Reference Guide

　あとは距離（距離行列で入力した場合は無視される）とクラスタリング方法を指定できます。選べるものの一覧はリファレンスを見てください。選択肢はいろいろ実装されているので、ここで「欲しいものがなくて困る」というシチュエーションはそう滅多にないと思います。

　linkageの返り値はnumpy配列です。細かいフォーマットについてはリファレンスを（ry。自分でこれをいじってどうこうしようという機会はあまりないと思うので、知らなくてもなんとなく使えます。リファレンスの説明ではだいたいZという変数に代入しているので、それに倣うと良いと思います。

fcluster

　クラスタリングしてくれます。

scipy.cluster.hierarchy.fcluster(Z, t, criterion='inconsistent', depth=2, R=None, monocrit=None)
scipy.cluster.hierarchy.fcluster — SciPy v1.3.0 Reference Guide

　いちいち細かく述べることはしませんが、Zにlinkageの返り値を入れ、tはスレッショルドでクラスタの分割/結合の基準として渡します。あとはcriterionに好きなアルゴリズムを選ぶ……という仕組みですね。

cophenet

　「cophenetic distances」なるものを計算します。

scipy.cluster.hierarchy.cophenet(Z, Y=None)
scipy.cluster.hierarchy.cophenet — SciPy v1.3.0 Reference Guide

　
　Y（クラスタリングの元になった距離行列。当然pdistフォーマット）を渡すとクラスタリング全体の評価指標を出してくれるので、その目的で使うのが良いと思います。返り値は2つあり得るんですが、Yを渡さなければ1つめ（全体の評価指標）は省略されて2つめだけ返ります。

　MATLABから輸入された関数らしく、MATLABのドキュメントを読んだ方がわかりやすいと思います。

コーフェン相関係数 - MATLAB cophenet - MathWorks 日本

　あとで書く実践編でちょっと触れます。

dendrogram

　デンドログラムを描きます。デンドログラム出しとけばとりあえず納得感があるので、これだけ出してみてからあれこれすれば良いと思います。

scipy.cluster.hierarchy.dendrogram(Z, p=30, truncate_mode=None, color_threshold=None, get_leaves=True, orientation='top', labels=None, count_sort=False, distance_sort=False, show_leaf_counts=True, no_plot=False, no_labels=False, leaf_font_size=None, leaf_rotation=None, leaf_label_func=None, show_contracted=False, link_color_func=None, ax=None, above_threshold_color='b')

scipy.cluster.hierarchy.dendrogram — SciPy v1.3.0 Reference Guide

　引数が泡吹いて倒れるほどいっぱいある。到底説明なんてできないので、代表的なものだけピックアップします。

　説明不要。linkageの結果を渡す。

　木を省略して表示するときのパラメータ。次のtruncate_modeと絡みます。

truncate_mode

　どのように木を省略するか。これはデータ数がある程度多いときに威力を発揮します。今回は使いません。

color_threshold

　色分けに絡むしきい値。木の中の最大ノード間距離との比率で色分けを決めます。

labels

　葉のラベル。

　matplotlibのAxesを渡せます。渡すとそこに描画されます。

　あとはまあ、いろいろ。基本的にはどれも表示フォーマットに絡む引数なので、見た目を変えたくなったら合う引数を探すという感じです。

実践編

　ではこれから実践してみます。

データを作る

　階層型クラスタリングは50>データ数くらいが適しています。多くても見づらいし計算量が嵩むからです。これくらいでのサイズの使いやすいデータはすぐに思い浮かびませんでしたが、sklearnのdigitsでラベルごとに平均を取ると良さげなことに気づいたので、そうします。

　ついでに可視化もする。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits

def gen_data():
    digits = load_digits()
    label_uniq = np.unique(digits.target)
    result = []
    for label in label_uniq:
        result.append(digits.data[digits.target == label].mean(axis=0))
    return result, label_uniq

def visualize():
    X, y = gen_data()
    fig, axes = plt.subplots(nrows=2, ncols=5)
    for ax, x, label in zip(axes.ravel(), X, y):
        ax.set_title(label)
        ax.imshow(x.reshape(8, 8))
    plt.savefig("data.png")

if __name__ == "__main__":
    visualize()

data.png

　よさげなので、この10件のデータでやってみます。

手法を選ぶ

　とりあえずmetric="euclidean"に固定してmethodを変化させ、評価指標を出してみましょう。

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, cophenet

def clustering_score():
    X, y = gen_data()
    methods = ["single", "complete", "average", "weighted",
               "centroid", "median", "ward"]
    for method in methods:
        S = pdist(X)
        Z = linkage(S, method=method)
        c, d = cophenet(Z, S)
        print("{0} {1:.3f}".format(method, c))

if __name__ == "__main__":
    clustering_score()

single 0.722
complete 0.752
average 0.769
weighted 0.766
centroid 0.681
median 0.730
ward 0.720

　average（UPGMA、いわゆる群平均法）がベストっぽいので、これを使うことにします。

　どれを選んでもだいたい0.75くらいのコーフェン相関係数なので、25%くらいはデータの性質が狂っていると思って結果を解釈する必要がある、ということだと思います。あまり細かいところを見ても無意味です。

クラスタに分ける

　fclusterを使うと適当な数のクラスタに分割できます。4つのクラスタに分割してみましょう。

from collections import defaultdict
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, fcluster

def clustering_fcluster():
    X, y = gen_data()
    S = pdist(X)
    Z = linkage(S, method="average")
    result = fcluster(Z, t=4, criterion="maxclust")
    d = defaultdict(list)
    for i, r in enumerate(result):
        d[r].append(i)
    for k, v in d.items():
        print(k, v)

if __name__ == "__main__":
    clustering_fcluster()

1 [1, 2, 3, 5, 8, 9] なんとなく形が似てるかな・・・？
2 [7]  # 独特の位置
3 [4, 6] # 似ているかと言うと微妙なものがある
4 [0]  # 独特の位置

　いまいち納得感の少ない結果になりました。

デンドログラムを描く

　デンドログラムを描いてみましょう。

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, dendrogram

def clustering_dendrogram():
    X, y = gen_data()
    S = pdist(X)
    Z = linkage(S, method="average")
    dendrogram(Z)
    plt.savefig("dendro1.png")

if __name__ == "__main__":
    clustering_dendrogram()

dendro1.png

　上の方でまとまっているのはあまりあてにならないので、実質的に意味があるのはとりあえず1,8と3,9あたりです。

　1,8は縦の真ん中あたりに集中するクラスタ。3,9は9の上の丸の左側を切り取ればほぼ3みたいな形になります。

　上のfclusterでやったのと同じ4つのクラスタに分かれるようにする方法がないのか？　基本的に違う考え方に基づいて着色されるので厳しいものがありますが、縦軸を見ながら4つに分かれるあたりを狙ってcolor_thresholdを決め打ちすると一応それに近いことはできます。

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, dendrogram

def clustering_dendrogram2():
    X, y = gen_data()
    S = pdist(X)
    Z = linkage(S, method="average")
    dendrogram(Z, color_threshold=31)
    plt.savefig("dendro2.png")

if __name__ == "__main__":
    clustering_dendrogram2()

dendro2.png

　このcolor_thresholdは距離の近さがスレッショルド未満のものを1つにまとめる、という挙動です。スレッショルドより大きい距離はabove_threshold_color（デフォルトは"b"で青）になります。今回はスレッショルドより上で1サンプルで別れてそのまま1クラスタを形成するという厄介な子がいるので微妙な結果になってしまいますが、もう少し性質の良いデータだとうまく合わせることはできると思います。

遊ぶ

　せっかくmethodが7種類あるので、それぞれでやってみます。

import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, cophenet, dendrogram

def clustering_dendro_each_method():
    X, y = gen_data()
    methods = ["single", "complete", "average", "weighted",
               "centroid", "median", "ward"]

    fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(10, 7))
    for ax, method in zip(axes.ravel(), methods):
        S = pdist(X)
        Z = linkage(S, method=method)
        c, d = cophenet(Z, S)
        dendrogram(Z, color_threshold=31, ax=ax)
        ax.set_title("{0} {1:.3f}".format(method, c))
    plt.savefig("dendro_each.png")

if __name__ == "__main__":
    clustering_dendro_each_method()

結果

　けっこう結果が変わるけど、細部はそれなりに揃っている感じ。コメントはとりあえず差し控えたいと思います。

まとめ

　このように使えます。便利関数や引数が多いのがややこしいだけで、基本的な使い方そのものは難しくないです。

　階層型クラスタリングは結果の良し悪しを客観的に評価するのが難しいので、割り切って使うと良いと思います。デンドログラムの細かい枝分かれの順序どうこうより、各クラスタの成因を判別分析や検定などでちゃんと調べてあげる方が本当は大切です。

scipy.interpolate.griddataの内挿方法による違いを比較

Fri, 14 Dec 2018 10:04:02 +0900

はじめに

　以前、3次元のサンプルデータを内挿してmatplotlibでうまくプロットする方法について記事にしました。

xyzの点データを内挿してmeshgridにしmatplotlibでプロットする - 静かなる名辞

　この記事では内挿のアルゴリズムをデフォルトのlinearにして使いましたが、他の方法ではどうなるのか気になったので実験してみました。

使えるアルゴリズム

　選択肢は3つだけです。

method : {‘linear’, ‘nearest’, ‘cubic’}, optional
Method of interpolation. One of
nearest
return the value at the data point closest to the point of interpolation. See NearestNDInterpolator for more details.
linear
tessellate the input point set to n-dimensional simplices, and interpolate linearly on each simplex. See LinearNDInterpolator for more details.
cubic (1-D)
return the value determined from a cubic spline.
cubic (2-D)
return the value determined from a piecewise cubic, continuously differentiable (C1), and approximately curvature-minimizing polynomial surface. See CloughTocher2DInterpolator for more details.
scipy.interpolate.griddata — SciPy v1.3.0 Reference Guide

　なんとなくcubicには1-Dと2-Dの2つがあって「1次キュービック補間と2次キュービック補間？　そんなのあったっけ」と思いがちですが、データが1次元か2次元かで使い分けられるだけで、ユーザが指定できるのは{‘linear’, ‘nearest’, ‘cubic’}のいずれかです。

　それぞれ

線形補間
最近傍補間
キュービック補間

　です。詳しい中身は知らなくても、いずれも名前くらいは聞いたことがあると思います。

実験

　二次元正規分布でサンプル数=128,512とし、それぞれの補間アルゴリズムで内挿します。結果をプロットして確認します。

　また、回帰とみなしてRMSEを出してみました。

　コードを以下に示します。

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats
from scipy import interpolate
from sklearn.metrics import mean_squared_error

def rmse(true, pred):
    return mean_squared_error(true.ravel(), pred.ravel())**(1/2)

def main():
    norm = stats.multivariate_normal(mean=[2.0, 3.0], cov=[[4, 2],[2,4]])

    # samples
    xy128 = np.random.uniform(low=-10, high=10, size=(128, 2))
    z128 = norm.pdf(xy128)
    xy512 = np.random.uniform(low=-10, high=10, size=(512, 2))
    z512 = norm.pdf(xy512)

    # xy meshgrid
    x = y = np.linspace(-10, 10, 500)
    X, Y = np.meshgrid(x, y)
    Z = norm.pdf(np.vstack([X.ravel(), Y.ravel()]).T).reshape(X.shape)

    # plot
    fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(10,5))
    plt.subplots_adjust(hspace=0.6, wspace=0.4)

    axes[0,0].pcolormesh(X, Y, Z, cmap="jet")
    axes[0,0].set_title("true data")
    axes[1,0].pcolormesh(X, Y, Z, cmap="jet")
    axes[1,0].set_title("true data")

    for i, (n_samples, xy, z) in enumerate(
            zip([128, 512], [xy128, xy512], [z128, z512])):
        axes[i,1].scatter(xy[:,0], xy[:,1], c=z, cmap="jet")
        axes[i,1].set_title("samples {}".format(n_samples))

        for j, i_method in enumerate(["nearest", "linear", "cubic"]):
            i_Z = interpolate.griddata(xy, z, (X, Y), method=i_method, 
                                       fill_value=0.0)
            axes[i,j+2].pcolormesh(X, Y, i_Z, cmap="jet")
            axes[i,j+2].set_title("{} {}\nrmse={:.5f}".format(
                i_method, str(n_samples), rmse(Z, i_Z)))
        
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　なお、RMSEを計算する都合上、fill_value=0.0としています。デフォルトはnanですが、それだと計算できないので……。一応実際にnanの状態でも確認し、nanになるのはグラフの端（このデータではほぼ0.0）の領域だけであることを確認して以上の判断をしました。

結果と考察

　プロットされた結果を示します。

result.png

　見ての通り、ダメダメな最近傍補間、まあまあな線形補間、群を抜いて良いキュービック補間という関係です。cubicにしておけば良いのでは？

　ただ、今回は真値そのものが補間が効きやすいなめらかなデータですが、実データはもう少しノイズが乗ったりして暴れることがあると思います。キュービック補間はオーバーシュート・アンダーシュートがあるらしいので、そういう場合でも対応できるように保険としてデフォルトがlinearになっているのかもしれません。まあ、無難なのはそっちでしょう。

　実用的には、両方やってみて大丈夫そうな方を選ぶことになるでしょう。

まとめ

　cubicがよかったです。

scipyで確率分布のサンプルと確率密度関数を生成する

Sun, 02 Dec 2018 03:56:08 +0900

　乱数データと確率密度関数を一緒にplotしてみたかったので、メモ。

概要

　scipy.statsでは様々な統計用のユーティリティが提供されています。大抵の分布はあるし、パラメータも好きに設定できます。

Statistical functions (scipy.stats) — SciPy v0.16.1 Reference Guide

　numpyにも充実したrandomモジュールがありますが、こちらは分布に従うデータの生成や、データのサンプリングなどしかできません（と思います）。

https://docs.scipy.org/doc/numpy/reference/routines.random.html

　なんとなく「データの生成はnumpyでできるけど、確率密度関数だとscipy使わないと駄目なのかな？」と思いがちですが、「データの生成」も実はscipyでできるので、numpyを使う必要性はありません。やったね。

方法

　この記事では正規分布でやってみます。

scipy.stats.norm — SciPy v0.16.1 Reference Guide

　scipy.stats.normというものを使うのですが、これは実はクラスに近い使い方ができます。

>>> from scipy import stats
>>> stats.norm
<scipy.stats._continuous_distns.norm_gen object at 0x7f8d839ede10>
>>> type(stats.norm)
<class 'scipy.stats._continuous_distns.norm_gen'>

　scipy.stats._continuous_distns.norm_gen objectというのはわかりづらいですが、早い話がファクトリであり、callするとオブジェクトを返します。

>>> stats.norm()
<scipy.stats._distn_infrastructure.rv_frozen object at 0x7f8da586ef28>

　素性はいまいちわかりませんが、rv_frozen objectという名前から想像できる通り、事実上「確率分布のインスタンス」のように使えます。

　つまり、分布のパラメータを渡してインスタンス化し、インスタンスを使ってそのパラメータの分布に従うサンプルを生成したり、確率密度関数を計算したりできるのですね。これを知ったときは喜びました。

　使えるメソッドはクラスによって微妙に違いがあるようですが、

rvs：Random Variates（確率変量）
pdf：Probability density function（確率密度関数）

　の2つが今回使いたい「分布に従うデータ」と「確率密度関数」を返すメソッドです。

　なお、これらのメソッドはstats.normのようなクラス（？）から呼ぶことも、インスタンス化したオブジェクト（？）（rv_frozenのオブジェクトと呼ぶべきか……）から呼ぶこともできますが、同じ名前でも使い方が違うので注意が必要です。

>>> from scipy import stats
>>> help(stats.norm.rvs)
Help on method rvs in module scipy.stats._distn_infrastructure:

rvs(*args, **kwds) method of scipy.stats._continuous_distns.norm_gen instance
    Random variates of given type.
    
    Parameters
    ----------
    arg1, arg2, arg3,... : array_like
        The shape parameter(s) for the distribution (see docstring of the
        instance object for more information).
    loc : array_like, optional
        Location parameter (default=0).
    scale : array_like, optional
        Scale parameter (default=1).
    size : int or tuple of ints, optional
        Defining number of random variates (default is 1).
    random_state : None or int or ``np.random.RandomState`` instance, optional
        If int or RandomState, use it for drawing the random variates.
        If None, rely on ``self.random_state``.
        Default is None.
    
    Returns
    -------
    rvs : ndarray or scalar
        Random variates of given `size`.

>>> help(stats.norm().rvs)
Help on method rvs in module scipy.stats._distn_infrastructure:

rvs(size=None, random_state=None) method of scipy.stats._distn_infrastructure.rv_frozen instance

　やりたいことはわかるし、こういう仕様にしたのも理解できるけど、私のような初心者は戸惑う。scipyのドキュメントの簡潔さも相まって。

実験

　こんなコードを動かしてみました。

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats

def main():
    normd = stats.norm(100, 7)
    x = np.arange(70, 130, 0.1)
    pdf = normd.pdf(x)
    samples = normd.rvs(1000)

    fig, ax1 = plt.subplots()
    ax1.hist(samples, bins=30, color="C0")
    ax2 = ax1.twinx()
    ax2.plot(x, pdf, color="C1")
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　2軸グラフの描き方についてはこちらを参考にしました。

[python]matplotlibで左右に２つの軸があるグラフを書く方法 - Qiita

　結果は、以下のようになりました。

result.png

　ちゃんとできていますね。

まとめ

　scipyでできることがわかった。scipyを信じていれば大抵のことはできる。

【python】sklearnのPCAでloading（主成分負荷量）を計算する

Sat, 31 Mar 2018 01:24:28 +0900

　PCA（主成分分析）のloading*1がほしいときがあります。

　sklearnでは一発では出ません。

　ドキュメントはここ。
sklearn.decomposition.PCA — scikit-learn 0.21.2 documentation

PCA.components_は確かにあるけど・・・
loadingを計算しよう
罠だった
共分散行列のときはどうするのか
loadingを使うと何が良いのか

PCA.components_は確かにあるけど・・・

　結論から先に言うと、PCA.components_はノルム1の固有ベクトルです。ノルムを測ってみましょう。

>>> import numpy as np
>>> from sklearn.datasets import load_iris
>>> from sklearn.decomposition import PCA
>>> iris = load_iris()
>>> pca = PCA(n_components=2)
>>> pca.fit(iris.data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.components_
array([[ 0.36158968, -0.08226889,  0.85657211,  0.35884393],
       [ 0.65653988,  0.72971237, -0.1757674 , -0.07470647]])
>>> np.linalg.norm(pca.components_, axis=1)
array([1., 1.])

　まあ、loadingも固有ベクトルには違いないのですが、ノルムを整えてやる必要があります。

loadingを計算しよう

　教科書などによく書いてあることですが、第主成分に対応する元の変数のloadingは次の式で出せます。

　は固有値。 eigenvectorは固有ベクトルで、元の変数の数だけ次元がありますから、これで良いわけです（雑な説明ですが・・・）。

　ということで、pythonで同様にやってみましょう。固有値はexplained_varianceに入っています。

>>> pca.components_*np.c_[np.sqrt(pca.explained_variance_)] # 縦ベクトルに変換する必要がある
array([[ 0.74322652, -0.16909891,  1.76063406,  0.73758279],
       [ 0.32313741,  0.35915163, -0.08650963, -0.03676921]])

　できました。これがloadingです。・・・と思ったけど、1を超えちゃってますね。なんてこった。

罠だった

　固有値は分散なので、データのスケールに依存します。

　とりあえず入力をスケーリングしてみよう。上の式は相関行列から行くときのものでした。なのでこれで平気なはず。

>>> from sklearn.preprocessing import StandardScaler as SS
>>> ss = SS()
>>> data = ss.fit_transform(iris.data)
>>> pca.fit(data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.components_*np.c_[np.sqrt(pca.explained_variance_)]
array([[ 0.89421016, -0.45081822,  0.99500666,  0.96822861],
       [ 0.35854928,  0.89132754,  0.02031465,  0.06299656]])

　1を超えなくなくてめでたし、ということよりも、数字が変わったことの方が問題で、これで本当に正しいのかという疑念が生じてきました。

　確認のために元の特徴と主成分の相関係数を直接測ってみます。

>>> X = pca.fit_transform(data)
>>> np.corrcoef(np.hstack([iris.data, X]), rowvar=False)
array([[ 1.00000000e+00, -1.09369250e-01,  8.71754157e-01,  8.17953633e-01,  8.91224479e-01,  3.57352114e-01],
       [-1.09369250e-01,  1.00000000e+00, -4.20516096e-01,  -3.56544090e-01, -4.49312976e-01,  8.88351481e-01],
       [ 8.71754157e-01, -4.20516096e-01,  1.00000000e+00,  9.62757097e-01,  9.91684422e-01,  2.02468206e-02],
       [ 8.17953633e-01, -3.56544090e-01,  9.62757097e-01,  1.00000000e+00,  9.64995787e-01,  6.27862218e-02],
       [ 8.91224479e-01, -4.49312976e-01,  9.91684422e-01,  9.64995787e-01,  1.00000000e+00,  2.08904471e-17],
       [ 3.57352114e-01,  8.88351481e-01,  2.02468206e-02,  6.27862218e-02,  2.08904471e-17,  1.00000000e+00]])

　下の二行の4列目までを見てください。微妙に誤差があるような気はしますが（小数点以下3桁でずれてきてるので微妙ってほど微妙でもないが）、たぶん同じ数字になっています。

　微妙な誤差については、丸め誤差などが累積した、実は計算間違ってる、といった理由が考えられます。前者ならまだ許せるけど、後者はやだな・・・。

共分散行列のときはどうするのか

　どうするんだろうね・・・。

　2019/06/22追記
　手順は増えますが、スケールを考慮すれば同様に行えるようです。

出典：http://manabukano.brilliant-future.net/document/text-PCA.pdf　p.10

loadingを使うと何が良いのか

　相関係数なので、「どれくらい効いてるか」がよくわかります。よくある0.3以下なら～とか0.7以上なら～という論法ができます。それだけといえばそれだけ。

　このように取扱が大変なので、固有ベクトルのまま解釈した方が楽かもという気がしてきました。各主成分の寄与率はexplained_variance_ratio_で得られる訳だし、寄与率の大きい軸の固有ベクトルの大きい次元を見ればどんな感じかはわかるし・・・。

　でもまあ、一応（入力をスケーリングすれば）大体出るということはわかったので、これでよしとします。

　共分散行列でやるときのやり方は、どなたか詳しい方に教えて頂けると助かります。

*1:主成分負荷量、あるいは因子負荷量とも（なぜか）言われますが、この記事ではloadingで通します。けっきょくヘタに和訳しないのがいちばんわかりやすい

【python】sklearnで因子分析を試す

Sat, 31 Mar 2018 00:22:11 +0900

　pythonで因子分析をやる人はあまりいないようだが、sklearnにはしっかりモデルが存在している。ついさっき気づいた。

sklearn.decomposition.FactorAnalysis — scikit-learn 0.20.1 documentation

　因子分析自体は前からどんなものなのか興味があり、かといってググるとRだったりSPSSだったりばっかり出てきて辟易していたのだが、sklearnにあると都合が良い。さっそく使ってみよう。

とりあえずirisをプロットする
とりあえずcomponentsを見る
使えることはわかった

とりあえずirisをプロットする

　私だけでも何十回もやってきた、世界中では何万回とやられてきたirisの二次元可視化をやってみる。

　次のようなコードを書いた。

# coding: UTF-8

from copy import deepcopy
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, FactorAnalysis as FA
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt

def decomp_and_plot(dataset, model, file_name):
    X = model.fit_transform(dataset.data)
    plt.figure()
    plt.scatter(X[:,0], X[:,1], c=dataset.target/len(dataset.target_names))
    plt.savefig(file_name)
    
def main():
    iris = load_iris()

    ss = StandardScaler()
    pca = PCA(n_components=2)
    pl = Pipeline([("scaler", ss), ("pca", deepcopy(pca))])
    fa = FA(n_components=2, max_iter=5000)

    decomp_and_plot(iris, pca, "pca_plt.png")
    decomp_and_plot(iris, pl, "spca_plt.png")
    decomp_and_plot(iris, fa, "fa_plt.png")

if __name__ == "__main__":
    main()

　PCA、変数をスケーリングしたPCA（相関行列を使うことと等価）、因子分析でそれぞれplotしてみる。

　結果はこれ。

PCAの結果

PCA(相関行列)の結果

　相関行列はぱっと見いまいち（この絵一枚でダメかどうかは判断できないが）。

因子分析の結果

　うーん、相関行列のとも違うし、なんとも言い難いというか、素人目にはぶっちゃけあんまり良くないように見えるのだが、確率モデルなのでノイズの存在を仮定して見るとこうなるということだろう。

とりあえずcomponentsを見る

　次のようなmain2を作り、実行した。

def main2():
    iris = load_iris()

    print(iris.feature_names)
    print("pca")
    pca = PCA(n_components=2)
    pca.fit(iris.data)
    print(pca.components_)

    print("fa")
    fa = FA(n_components=2, max_iter=5000)
    fa.fit(iris.data)
    print(fa.components_)

　結果

['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
pca
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [ 0.65653988  0.72971237 -0.1757674  -0.07470647]]
fa
[[ 0.72577591 -0.17754023  1.75733754  0.73196365]
 [-0.37036948 -0.24060118  0.02793388  0.04121372]]

　プロット結果から予想される通り、両者のcomponentsはよく似通っている。

　これがloadingなのかどうかはぶっちゃけよくわからないのだが（というか1を超えてくる時点でたぶん違うのだろうが）、とりあえずloadingだと思って解釈する。

　第一因子は花弁の長さと幅、がく片の長さに対応しているので花の大きさに対応しているっぽい。花の大きさとがく片の幅はなぜか若干反比例する。

　第二因子は花弁に関する係数が小さいので、がく片の大きさを表す因子と言って良さそうである。

　こんなところか。

使えることはわかった

　だから何？　って言われると、正直答えに窮しますが・・・とにかく使えます。主成分分析で良いじゃんと言われたら何も言い返せません。
　

【python】pythonで主成分分析のバイプロット

Wed, 28 Mar 2018 23:13:05 +0900

　バイプロット（Biplot）という主成分分析（PCA）の結果の可視化方法があります。

　すごく大雑把に言うと、PCAによる写像の前の空間の各特徴（軸）が写像先の空間のどこに向いているかを可視化する方法です。

　具体的には、主成分ベクトル（因子負荷量などを使う場合もあります）と散布図を同じ図にplotします。これらを組み合わせることで、元の空間の性質が二次元（もしかしたら3次元）で手に取るようにわかります*1。

　バイプロットはR言語だと簡単に描けるらしいのですが、我らがpythonには（少なくとも一般的なライブラリには）そんな便利なものはありません。ちょっと困るのですが、幸い英語圏にはちらほらやりかたの情報があります。しかし、それはそれでページごとにやってることが違ったりして、（申し訳ないのですが）微妙に信用できなかったりします。

　で、けっきょく自分で書いてみることにしました。なお、参考にしたのはこの辺です。

方針

　まずsklearnの公式ドキュメントをできるだけ良く読み込みます。

sklearn.decomposition.PCA — scikit-learn 0.22.1 documentation

　PCA.components_が固有ベクトルであり、データをセンタリングしてこれと掛けるとPCAの出力が出てくることは前回の記事で確認しました。

　固有ベクトル行列が主成分*元のデータの特徴という形になっているとして、横に見ると負荷量（みたいなもの。本当は対応する固有値のsqrtを掛け算してやらないといけない）に、縦に見ると元の写像先で表現された特徴の軸になります。

　つまり、その軸をプロットするだけです。

　なお、この辺は微妙に議論があるようです。私もこれがどこまで正しい方法なのかは自信が持てません。

　参考：
色々と考えてみる: 文系のための「主成分分析の可視化」（２）

　だけど今回は、データをセンタリングしてPCAを学習させた上で、各軸に対応するone-hot vectorを渡してtransformしたら確かに上に書いた方法通りで上手く行きました（biplotの線の上に載った）。なので、「これで良いんだろう」と勝手に判断しました。どこまで妥当かはよくわからないんですけど。

実装

　こんな感じで書きました。

# coding: UTF-8

from sklearn.datasets import load_iris, load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

def biplot(dataset, scale=False, arrow_mul=1, text_mul=1.1):
    if scale:
        ss = StandardScaler()
        X = ss.fit_transform(dataset.data)
    else:
        X = dataset.data

    if hasattr(dataset, "feature_names"):
        feature_names = list(dataset.feature_names)
    else:
        feature_names = ["F{0}".format(i)
                         for i in range(dataset.data.shape[1])]

    pca = PCA(n_components=2)
    X = pca.fit_transform(X)

    x_data = X[:,0]
    y_data = X[:,1]

    pc0 = pca.components_[0]
    pc1 = pca.components_[1]

    plt.figure()
    plt.scatter(x_data, y_data,
                c=dataset.target/len(set(dataset.target)),
                marker=".")

    for i in range(pc0.shape[0]):
        plt.arrow(0, 0, 
                  pc0[i]*arrow_mul, pc1[i]*arrow_mul,
                  color='r')
        plt.text(pc0[i]*arrow_mul*text_mul,
                 pc1[i]*arrow_mul*text_mul,
                 feature_names[i],
                 color='r')
    plt.show()

def main():
    iris = load_iris()
    wine = load_wine()

    biplot(iris, arrow_mul=2.5, scale=True)
    biplot(wine, arrow_mul=6, scale=True)

if __name__ == "__main__":
    main()

　今回はsklearnのデータセットを渡す形で関数にまとめました。ま、もしこのコードを流用したい人がいたら、必要なロジックだけ上手く切り出してください。

　結果は、こんな画像が出ます。

irisのバイプロット

wineのバイプロット

　上手く行ってる感じです。

　なお、上のコードでは変数をスケーリングしています（相関行列でPCAするのと等価）。スケーリングしなくてもできますが、やった方が矢印の長さが揃いやすいです（逆に変数のスケールを重視してPCAしたいときは、スケーリングしてはいけない。ケースバイケース）。

まとめ

　これくらい自作しなくても済めば良いのにと思いました。

*1:本当に手に取るようにわかるかはデータと見る人に依存しますが・・・

【python】numpyで主成分分析を実装してみた

Wed, 28 Mar 2018 22:21:01 +0900

　numpyでPCA（principal component analysis：主成分分析）を実装してみました。自分の理解を深めるためです。

　sklearnに実装されているものと同じ結果を出すことを目標にしました。最終的には上手く行きました。

概要
実装
結果
まとめ

概要

　主成分分析のアルゴリズムの解説は他に譲ります。これは実装してみた記事です。

　実装のやり方は色々あるようですが、一番基本的な（だと思う）共分散行列の固有値と固有ベクトルを求める方法で行きます。

　やるべきこととしては、

データをセンタリングする（列ごとに平均を引く）
共分散行列を計算する
固有値と固有ベクトルを計算
データを固有ベクトルを使って写像する

　これらを実装すれば行けるはずです。というか、これで行くことにしました。

実装

　書いたソースコードを以下に示します。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

class MyPCA:
    def __init__(self, n_components=2):
        self.n_components = n_components

    def fit_transform(self, X):
        """横着してfit_transformしか実装してない
        """

        # 平均を0にする
        X = X - X.mean(axis=0)

        # 共分散行列を作る
        self.cov_ = np.cov(X, rowvar=False)
        
        # 固有値と固有ベクトルを求めて固有値の大きい順にソート
        l, v = np.linalg.eig(self.cov_)
        l_index = np.argsort(l)[::-1]
        self.l_ = l[l_index]
        self.v_ = v[:,l_index] # 列ベクトルなのに注意

        # components_（固有ベクトル行列を途中まで取り出す）を作る
        self.components_ = self.v_[:,:self.n_components].T

        # データとcomponents_をかける
        # 上と下で二回転置してるのアホ・・・
        T = (np.mat(X)*(np.mat(self.components_.T))).A

        # 出力
        return T

def main():
    iris = load_iris()

    pca = PCA(n_components=2)
    sklearn_X = pca.fit_transform(iris.data)

    my_pca = MyPCA()
    my_X = my_pca.fit_transform(iris.data)

    print(pca.explained_variance_)
    print(my_pca.l_)

    print(pca.components_)
    print(my_pca.components_)

    plt.figure()
    plt.scatter(sklearn_X[:,0], sklearn_X[:,1], c=iris.target/3)
    plt.savefig("sklearn_resut.png")

    plt.figure()
    plt.scatter(my_X[:,0], my_X[:,1]*-1, c=iris.target/3)
    plt.savefig("my_result.png")

if __name__ == "__main__":
    main()

　numpyを使ったので簡単に書けました。アルゴリズム部分はコメントで解説を入れたので、それを読めばどんな感じかは理解して頂けると思います。

結果

　mainのテキスト出力を見ると、次のようになっていました。

# 固有値
[4.22484077 0.24224357]
[4.22484077 0.24224357 0.07852391 0.02368303]

# components_
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [ 0.65653988  0.72971237 -0.1757674  -0.07470647]]
[[ 0.36158968 -0.08226889  0.85657211  0.35884393]
 [-0.65653988 -0.72971237  0.1757674   0.07470647]]

　固有値が余計に出ちゃってますが、これは別に構いません。また、componentsの2次元目が符号反転していますが、これも特に問題ないこと（のはず）なので無視します。

　自作の方は第二主成分を反転させてプロットしてみました。

sklearnのPCAでirisを可視化

自作PCAでirisを可視化

　同じ図を2つ載せるなって怒られそうですが・・・とにかく上手く行ったようです。

まとめ

　numpyで実装してみたら思ったより簡単だったので、これで当分は「わかった気」になれそうです。

　ただ、今回は特異値分解やらなかったので、それはまた宿題ということで・・・。

【python】scipy.statsのzscoreで警告が出るときの対策

Tue, 20 Mar 2018 17:30:20 +0900

概要

　z得点を計算しようとしたとき、このような警告を見かけることがあります。

RuntimeWarning: invalid value encountered in true_divide

　これが出た場合、結果にはnanが含まれています。なので後段の分析で落ちたりします。

>>> import numpy as np
>>> from scipy.stats import zscore
>>> a = np.array([[1,2,3,4], [1,2,3,4],[1,3,4,5]])
>>> zscore(a, axis=0)
stats.py:2248: RuntimeWarning: invalid value encountered in true_divide
  return (a - mns) / sstd
array([[        nan, -0.70710678, -0.70710678, -0.70710678],
       [        nan, -0.70710678, -0.70710678, -0.70710678],
       [        nan,  1.41421356,  1.41421356,  1.41421356]])

　どうしてエラーになるかというと、z得点は標準偏差で割るので、標準偏差が0だと0除算エラーが発生するからです。標準偏差0の列が含まれるようなゴミだらけの汚い疎行列をそのまま入れると、これが出来ます。

対策

　どうせ標準偏差0の軸とか要らないので、予め消し飛ばしておく。

>>> a[:,np.std(a, axis=0) != 0]
array([[2, 3, 4],
       [2, 3, 4],
       [3, 4, 5]])
>>> zscore(a[:,np.std(a, axis=0) != 0], axis=0)
array([[-0.70710678, -0.70710678, -0.70710678],
       [-0.70710678, -0.70710678, -0.70710678],
       [ 1.41421356,  1.41421356,  1.41421356]])

　めでたしめでたし。

【python】混合ガウスモデル (GMM)でハード・ソフトクラスタリング

Tue, 06 Mar 2018 04:39:50 +0900

はじめに

　先日はFuzzy c-meansによるソフトクラスタリングを行いました。

【python】skfuzzyのFuzzy c-meansでソフトクラスタリング - 静かなる名辞

　ソフトクラスタリングの有名な手法としてはFuzzy c-meansの他に、混合ガウスモデル（混合正規分布モデル）を使った手法があります。この手法はデータが「複数の正規分布から構成されている」と仮定し、その正規分布のパラメタ*1をEMアルゴリズム（expectation–maximization algorithm）という手法を使って最尤推定します。

　ごちゃごちゃと書きましたが、要するに「3つのクラスタにクラスタリングしたければ、（各クラスタのデータの分布が正規分布に従うと仮定して）3つの正規分布が重なりあってると思ってGMMを使って解く」という乱暴なお話です。正規分布が重なりあっているとみなすということは、どの分布に属するかも確率でわかる訳で、これがソフトクラスタリングに使える理由です。ハードクラスタリングに使いたいときは、確率最大のクラスタラベルに振ることになるかと思います。

　このGMM、pythonではsklearnに入っているので簡単に使えます。

sklearn.mixture.GaussianMixture — scikit-learn 0.20.1 documentation

　ということで、他のクラスタリング手法と比較してみることにしました。

実験の説明

　先日の記事でやったのと同様、irisをPCAで二次元に落としたデータに対してクラスタリングを行います。クラスタリング結果（所属するクラスタの確率）はirisが3クラスのデータなのを利用し、色（RGB）で表現します。

　比較するクラスタリング手法はk-means（ハード）、Fuzzy c-means（ソフト）、GMM（ハード・ソフト）です。

　前回はFuzzy c-meansのパラメタmを動かして結果を見たりしましたが、今回これは2で決め打ちにします。

　実験用ソースコードは次のものです。走らせるにはいつもの定番ライブラリ以外にscikit-fuzzyというライブラリを入れる必要があります（あるいはFuzzy c-means関連の部分をコメントアウトするか。でもskfuzzyはpipで一発で入るし、入れておいても別に損はない）。

# coding: UTF-8

import numpy as np

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans as KM
from sklearn.mixture import GaussianMixture as GMM
from matplotlib import pyplot as plt

from skfuzzy.cluster import cmeans

def target_to_color(target):
    if type(target) == np.ndarray:
        return (target[0], target[1], target[2])
    else:
        return "rgb"[target]

def plot_data(data, target, filename="fig.png"):
    plt.figure()
    plt.scatter(data[:,0], data[:,1], c=[target_to_color(t) for t in target])
    plt.savefig(filename)

def gen_data():
    iris = load_iris()
    pca = PCA(n_components=2)
    return pca.fit_transform(iris.data), iris.target

def main():
    data, target = gen_data()
    plot_data(data, target, filename="origin.png")

    km = KM(n_clusters=3)
    km_target = km.fit_predict(data)
    plot_data(data, km_target, filename="kmeans.png")

    cm_result = cmeans(data.T, 3, 2, 0.003, 10000)
    plot_data(data, cm_result[1].T, filename="cmeans_2.png")

    gmm = GMM(n_components=3, max_iter=1000)
    gmm.fit(data)
    gmm_target = gmm.predict(data)
    gmm_target_proba = gmm.predict_proba(data)
    plot_data(data, gmm_target, filename="gmm.png")
    plot_data(data, gmm_target_proba, filename="gmm_proba.png")

if __name__ == "__main__":
    main()

結果

オリジナルデータ

元データ

　これが元のデータです。できるだけこれに近いようなクラスタリング結果を得ることを目標とします。

k-means

　図の左側のクラスタは分離できていますが、右側は割と悲惨です。クラスタ同士が隣接していて細長い形だったりすると上手く行かないことが多いのがk-meansの特徴です。

c-means

Fuzzy c-means

　こうして見るとc-meansは「ファジー理論を入れて境界を曖昧にしたk-means」という気がしてきます。実際アルゴリズムもそんな感じなんですけど。

GMM

GMM-based clustering (hard)

GMM-based clustering (soft)

　一見して「おお」って感じですね。k-means、c-meansと比較して、元データのラベルに近いクラスタリング結果が得られています（図の右側の2つのクラスタの境界が右肩上がりになっている）。まあ、ちょっと元データのラベルとはずれているんですが（右下の方はかなり怪しい）、普通はこちらの方がk-meansやc-meansより「良い」クラスタリング結果だ、と判断されることが多いでしょう。

　どうしてこうなるのかというと、「irisのデータが正規分布していた」ということに尽きます。ま、アヤメの花びらの大きさとかのデータですから、正規分布しているんでしょう、きっと。

　こうして見るとGMMの方が良さそうな気もしますが、「ちゃんと正規分布してるか」が怪しいとちょっと適用するのを躊躇うのと、あと計算コスト自体はk-meansより高いはずなので*2、いまいちk-meansと比べて使われていない、というのが実情に近いかもしれません。

まとめ

　GMMを使ってみたらけっこう良かったです。

*1:一次元なら平均と分散、多次元なら共分散みたいな話になってくるのだろうか？

*2:Fuzzy c-meansとどっちが良いかは未調査

【python】正準相関分析（Canonical Correlation Analysis）を試してみる

Fri, 16 Feb 2018 02:13:14 +0900

　正準相関分析を使うと、2つの多次元データ同士の関連性を分析できるらしい。

　面白そうなので試してみた。ちなみに正準相関はsklearn.cross_decomposition.CCAで使える。正準相関自体の解説はほとんどしないので、文中のリンクを参考にして欲しい*1。

一応概要だけ
ノイズに埋もれた波形を取り出す
もうちょっとデータ分析っぽいことをしてみる

一応概要だけ

　代表的な多変量解析の手法（といって良いのかどうか少し悩むけど）として、主成分分析や重回帰分析が存在する。

主成分分析：一つの多変量データを直交するより少ない変数に縮約する
重回帰分析：一つの多変量データを一つの単変量データに変換する

　主成分分析にしろ重回帰分析にしろ、変換の係数だったり行列だったりを求めてそれで変換するのが実際にやることである。

　さて、正準相関は上の流れで説明すると、

正準相関分析：二つの多変量データをそれぞれ直交するより少ない変数に縮約して、かつ二つの変換されたデータの間で相関を最大化する

　という目的の分析である。主成分分析と重回帰を混ぜた感じ。

　気づいた人もいると思うけど、多変量vs多変量のデータでどちらかを単変量に分解して個別に重回帰で解くことも可能である。それに対するメリットとしては、

個別に重回帰するより全体の構造みたいなものを捉えられる可能性がある
個別に重回帰すると係数の数が全体でとても多くなるので解釈が面倒くさいが、一度次元を下げて直交した空間に持っていくことでそこが楽になる

　というあたりがあり、要するに解釈性がいいということ。

　この説明でもよくわからん、という人は、ニューラルネットのオートエンコーダーとか思い浮かべていただくとかえってわかりやすいかもしれない。

ノイズに埋もれた波形を取り出す

　参考URLの通りにやることにする。

　単一の信号源に複数のプローブを当てていて、それぞれに独立のノイズが乗って信号が埋もれてしまった・・・みたいな状況から元の信号を取り出そうとしているらしい。脳波計測とかで使えるのだろうか？

　参考URL：https://www.jstage.jst.go.jp/article/jnns/20/2/20_62/_pdf

　とりあえずこのようなコードを書き、

# coding: UTF-8
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cross_decomposition import CCA

def plot_wave(data, filename):
    fig, [ax1,ax2] = plt.subplots(2,1,figsize=(16,9))

    ax1.plot(data[:,0], color="b")
    ax2.plot(data[:,1], color="r")

    plt.savefig(filename)

def gen_pulse_data():
    common_pulse = np.array([-1]*50 + [0]*50 + [1]*50 + [0]*50 + 
                            [-1]*50 + [0]*50 + [1]*50 + [0]*50, dtype=np.float64)
    common_pulse += (np.random.random(common_pulse.shape) - 0.5)*0.1

    noise1 = (np.random.random(common_pulse.shape) - 0.5)*50
    noise2 = (np.random.random(common_pulse.shape) - 0.5)*50

    data1 = np.vstack([common_pulse + noise1,  common_pulse - noise1]).T
    data2 = np.vstack([common_pulse + noise2,  common_pulse - noise2]).T

    return data1, data2

def main():
    X1, X2 = gen_pulse_data()
    plot_wave(X1, "X1.png")
    plot_wave(X2, "X2.png")

    cca = CCA(n_components=2)
    cca.fit(X1, X2)

    Y1 = cca.transform(X1)
    Y2 = cca.transform(X2)
    plot_wave(Y1, "Y1.png")
    plot_wave(Y2, "Y2.png")

if __name__ == "__main__":
    main()

　実行する。

　まずは元の信号。

　わかる訳ねえな、という感じ。

　CCAでX1とX2の相関が最大になるような変換を計算し、その変換に基いてX1を変換したものをY1とすると、

　こんな感じになった。Y2も同じようなものなので省略。ここでは2次元出しているが、元のパルス信号が1次元なので2次元目（下）はノイズだけ出ている。

　まあ、上手く動いているのではないだろうか。

もうちょっとデータ分析っぽいことをしてみる

　如何せん、「ノイズに埋もれた信号を取り出せる！」というだけでは、データ分析っぽくなくて（個人的には）面白くない。正準相関自体はもっと色々なことができるはず。

　ここで足を引っ張るのは「正準相関向きのサンプルデータが見つからない」ということである。

正準相関向きのデータを探すのは困難

　2つの多次元データが対応しているようなデータで、適当にわかりやすいものがあれば良いのだが・・・なかなか良いデータがない。上に挙げた解説論文でも、「知名度は低い」とか書かれちゃってるし、正準相関自体、ニッチな感じがする。そこが素敵なのだが。

　一応、ネット上にある解説例だと、

統計学入門−第19章

　医学の分野で、肝機能の検査値（複数）と腎機能の検査値（複数）の対応を見るとか、

http://ogasun.la.coocan.jp/hanbetsubunseki.pdf

　中学生の体格（身長、体重、座高とか）と運動能力（50m走、走り幅跳びとか）の対応を見るとか、

　そういう感じのことをやっているのだが、この手のデータを探してくるのがまず面倒くさいし、見つけてもプログラムに流し込めるようにするまでがまた苦行だろうな、ということは容易に想像できるのである。

　この点で悩んで、この記事も一週間くらい出すか出さないか迷ってたんだけど、やることにした。ただ、結局良いデータは見つからなかったので、それっぽくでっちあげることにした。

作成したデータ

　ある架空の中学校で集計したという設定の、20人の生徒のデータである。「学外での勉強や取り組み」と「学校の成績」が対応付いている。

　「学外での勉強や取り組み」には、

一ヶ月に何冊読書するか
一年に何回博物館に行くか
毎週何日塾に通っているか
毎日何時間自習しているか

　の4つの変数がある。一方、「学校の成績」は、

国語
数学
社会
理科
英語

　の5つの科目があり、5段階評価で成績が付く。

　本来であれば適当に線形モデルでも作ってあげて数字を作るべきところだが、面倒くさいので私の想像で適当に埋めた（ツッコミポイント）。

　一応、次のような方針を考え、それに沿ったデータになるようにでっちあげた・・・つもり。

読書量と国語の成績は比例する
博物館に行った回数と社会、理科の成績は比例する
塾に通う頻度、自習時間は成績全体に影響を及ぼす

　よって、こういう結果が出てくるか、という勝負になる。

実験と結果

　こういうプログラムを書いた。

# coding: UTF-8
import numpy as np
from scipy.stats import pearsonr
from sklearn.cross_decomposition import CCA

def gen_data():
    # X1:
    # 毎月何冊の本を読むか,
    # 一年に何回博物館に行くか,
    # 塾に週何日通うか,
    # 毎日何時間自習するか
    X1 = np.array([[1,0,2,1],
                   [3,2,4,2],
                   [0,0,2,0],
                   [9,4,2,1],
                   [1,1,3,1],
                   [8,1,6,3],
                   [0,9,7,8],
                   [2,2,4,1],
                   [5,0,0,1],
                   [2,0,4,0],
                   [0,0,7,8],
                   [4,4,2,2],
                   [5,1,2,1],
                   [1,1,5,2],
                   [8,6,2,1],
                   [0,0,0,1],
                   [6,1,3,1],
                   [2,0,3,1],
                   [4,8,5,3],
                   [5,0,1,1]])

    # X2:
    # 国語,数学,社会,理科,英語の成績
    X2 = np.array([[3,3,3,3,3],
                   [4,3,4,4,5],
                   [2,2,3,3,2],
                   [5,4,3,3,3],
                   [3,3,4,4,4],
                   [5,5,5,4,5],
                   [3,5,5,4,5],
                   [4,4,4,5,3],
                   [5,3,3,3,3],
                   [3,4,3,4,3],
                   [5,5,4,5,5],
                   [4,4,5,5,3],
                   [4,3,3,3,3],
                   [4,4,5,4,5],
                   [5,3,5,5,3],
                   [2,2,2,1,2],
                   [5,3,4,4,4],
                   [3,4,3,4,3],
                   [5,5,5,5,5],
                   [5,3,3,3,3]])
    return X1, X2

def main():
    X1, X2 = gen_data()

    cca = CCA(n_components=4)
    cca.fit(X1, X2)

    print("Correlation Coefficient")
    for i in range(4):
        print("{0}:{1:.4f}".format(i, pearsonr(cca.x_scores_[:,i], cca.y_scores_[:,i])[0]))

    print("")
    np.set_printoptions(formatter={'float': '{: 0.4f}'.format})
    print("X1 loadings")
    print(cca.x_loadings_.T)
    print("")
    print("X2 loadings")
    print(cca.y_loadings_.T)

if __name__ == "__main__":
    main()

　「学外での勉強や取り組み」=X1と「学校の成績」=X2を4次元の空間上に写像して相関を最大化する、という問題を解かせる。軸同士は直交していて無相関なので、写像したデータの各軸の値同士の相関だけ見てやれば良い。写像したデータは、cca.x(or y)_scores_かcca.transform(X1(or X2))で取得できる*2。

　あとはX1とX2の各成分が、写像先の各軸にどれだけ寄与しているかがわかれば良い。そのためにはcca.x(or y)_loadings_を見る。転置した方が見やすいのでそうしている。

　こうして見ると、PCAに似ている。実際、CCAとPCAは親戚らしい。ま、あまり理論的な話に深入りしてもボロが出るので、これくらいにしておく。

　さて、結果はこのようになった。

Correlation Coefficient
0:0.9558
1:0.8978
2:0.5980
3:0.2927

X1 loadings
[[-0.4224  0.4244  1.0047  0.7353]
 [ 0.9326  0.4315  0.2511  0.3450]
 [ 0.2350  0.9269 -0.1714 -0.2514]
 [-0.4357  0.4062 -0.0634  0.8007]]

X2 loadings
[[-0.0558  0.7802  0.6574  0.6206  0.8044]
 [ 0.9392  0.5254  0.4702  0.3366  0.4271]
 [-0.2184 -0.0705  0.9353  0.5442 -0.3528]
 [-0.4310  0.0676 -0.2136 -0.8751  0.0008]]

　まず見るべきはCorrelation Coefficientで、写像先の空間の軸にどれだけ相関（＝やった意味）があるかを示している。0,1,2次元目はまあまあ強い相関だが、3次元目は相関係数0.3じゃ大した意味はなさそうだな、という風に解釈しておく。

　次にX1 loadingsとX2 loadingsを見る。X1 loadingsは4*4、X2 loadingsは4*5で、つまり行が写像先の軸、列が元の空間の軸に対応するように表示している。

　X1 loadingsの各行を見ていくと、

1行目

　塾と自習に熱心

2行目

　読書

3行目

　博物館

4行目

　自習と博物館だけ？

　なんとか解釈できる。数字がでかいところだけ重視するのがこつ。X2 loadingsも同様にやると、

1行目

　国語以外のすべて。国語にはほぼ中立。特に強いのは英語

2行目

　国語。他もそれなりに

3行目

　社会と理科

4行目

　理科にとてもネガティブ。全体的にネガティブな感じ

　ここまで出揃えば後はなんとかなる。このデータを作った方針を再掲する。

読書量と国語の成績は比例する
博物館に行った回数と社会、理科の成績は比例する
塾に通う頻度、自習時間は成績全体に影響を及ぼす

　0次元目は「塾に通う頻度、自習時間は成績全体に影響を及ぼす」とに、1次元目は「読書量と国語の成績は比例する」に、2次元目は「博物館に行った回数と社会、理科の成績は比例する」に対応していることがわかり、まあ妥当な結果なんじゃないの、という気はする。相関係数の低い3次元目はそこまで重視する必要はない。

　今回は先に方針を決めてデータをでっち上げたのであまり感動がないような気もするが、実際はデータにどんな構造があるのかは分析してみないとわからない。その構造を理解する上で正準相関が役に立つことは、上の例でなんとなく理解できた。

*1:正準相関でググって1ページ目に出てくるようなページばかり・・・

*2:今回はどちらも同じ値が返るが、transformだと学習時とは違うデータも入れられる

【python】95%信頼楕円/確率楕円を描画する

Wed, 14 Feb 2018 23:55:00 +0900

　「ライブラリあるやろｗ」と思ったら、なくて顔面蒼白になった。

　しょうがないから調べて実装した。

理論的なもの

　ちゃんと数式を書いて説明する気概がないので、アバウトに説明する。

　適当な二次元正規分布のデータがあるとする。PCAと同じ要領で分散共分散行列を対角化する。

　対角化した行列の対角成分（＝固有値）は、データを軸同士の相関がないような空間に変換して（＝要するにぐるっと回して）あげたときの変換先の軸上における分散である。

　分散がわかれば、一次元のときの信頼区間的なものがわかる。それを決めるためのデータの中心位置からの距離にはマハラノビス距離という名前が付いている。そして、これのニ乗は分布になるので、けっきょく累積確率だけ決めれば適当に定まる、ということがわかる*1。相関がなくなるように軸を作ったので、各軸ごとで求めてやれば（その軸上での）信頼楕円の径がわかる。

　この段階で楕円の幅と高さがわかっていることになるので、あとは適当に回転角を計算すると、

楕円の中心位置（単純に全データの平均で良い）
幅と高さ
回転角

　
　がわかることになり、楕円が描ける。

実装

　実装したものをそのまま貼っておきます。説明はしないので、参考にしたければしてください。

　一応入力データを色々変えてみておかしくなさそうなことは確認していますが、ミスがないとは言い切れません。ご注意ください。見つけたら指摘して頂けると嬉しいです。

# coding: UTF-8

import numpy as np
from scipy.stats import chi2

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

class ConfidenceEllipse:
    def __init__(self, data, p=0.95):
        self.data = data
        self.p = p

        self.means = np.mean(data, axis=0)
        self.cov = np.cov(data[:,0], data[:,1])

        lambdas, vecs = np.linalg.eigh(self.cov)
        order = lambdas.argsort()[::-1]
        lambdas, vecs = lambdas[order], vecs[:,order]

        c = np.sqrt(chi2.ppf(self.p, 2))
        self.w, self.h = 2 * c * np.sqrt(lambdas)
        self.theta = np.degrees(np.arctan(
            ((lambdas[0] - lambdas[1])/self.cov[0,1])))
        
    def get_params(self):
        return self.means, self.w, self.h, self.theta

    def get_patch(self, line_color="black", face_color="none", alpha=0):
        el = Ellipse(xy=self.means,
                     width=self.w, height=self.h,
                     angle=self.theta, color=line_color, alpha=alpha)
        el.set_facecolor(face_color)
        return el

def gen_data():
    return np.random.multivariate_normal([3,3], [[0.3,-0.2],[-0.2,1]], size=100)
def main():
    data = gen_data()

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.scatter(data[:,0], data[:,1], color="b", marker=".", s=3)

    el = ConfidenceEllipse(data, p=0.95)
    ax.add_artist(el.get_patch(face_color="blue", alpha=0.5))
    plt.savefig("img.png")

if __name__ == "__main__":
    main()

　こんな感じの絵が出てきます。

　正しいはずだけど、詳しくチェックはしていないので、自己責任でご利用ください。

参考文献

*1:てきとーすぎる説明だな

統計 - 静かなる名辞

【python】matplotlibのboxplotで外れ値を表示しないようにする

はじめに

シンプルに表示させない

外れ値の計算そのものをやめる

まとめ

【python】相関係数行列をstatsmodelsを使って描く

はじめに

使い方

まとめ

pythonで相関係数を計算する方法いろいろ3種類

はじめに

データの確認

numpyでやる

pandasでやる

scipyを使う

あと思ったこととか

まとめ

scikit-learnで目的変数を対数変換したりするTransformedTargetRegressor

はじめに

目的変数の対数変換

やってみる

おまけ：説明変数を対数変換したいとき

まとめ

ロジスティック回帰が線形分離不可能な分類問題を解けないことの説明

はじめに

y=0.5を代入すればいい

非線形の問題も解く方法

コサイン距離は距離じゃないんだから、勘違いしないでよねっ！

き、記事タイトルに意味なんてないんだからねっ！

コサイン類似度のことくらい自分で調べなさいっ！

距離の定義を知らないの？ しょ、しょうがないから教えてあげるわ

わからないの？ ……ばか

距離として扱うと困るのかって？ ……困るに決まってるじゃないっ、わからずや！

「じゃあどう呼べば良いのか」って？ そんなの自分で考えてよね！

わかったなら感謝しなさい。……え、ありがとう？ べ、べつに喜ばれても嬉しくなかんないんだからっ！

【python】ロジスティック回帰で確率値で学習させる

はじめに

データ

scikit-learnでは（たぶん）できない

statsmodelsでやる

まとめ

【python】statsmodelsでt検定する方法

はじめに

使い方の概要

実際にやってみる

まとめ

【python】scipyでt検定する方法まとめ

概要

t検定の概要

使う関数

引数と注意点

やってみる

まとめ

statsmodelsでやる方法

【python】scipyで階層型クラスタリングするときの知見まとめ

はじめに

関数がいっぱいある

使い方

linkage

fcluster

cophenet

dendrogram

実践編

データを作る

手法を選ぶ

クラスタに分ける

デンドログラムを描く

遊ぶ

まとめ

scipy.interpolate.griddataの内挿方法による違いを比較

はじめに

使えるアルゴリズム

実験

結果と考察

まとめ

scipyで確率分布のサンプルと確率密度関数を生成する

概要

方法

実験

距離の定義を知らないの？　しょ、しょうがないから教えてあげるわ

わからないの？　……ばか

距離として扱うと困るのかって？　……困るに決まってるじゃないっ、わからずや！

「じゃあどう呼べば良いのか」って？　そんなの自分で考えてよね！

わかったなら感謝しなさい。……え、ありがとう？　べ、べつに喜ばれても嬉しくなかんないんだからっ！