VBAチートシート

これは自分用のチートシートで随時更新型エントリです。

データの扱いはSQLやPythonですることが多いですが、やっぱり書類や設計書作成ではVBAは強い味方です。

使い方

Alt+F11でVBAを起動
標準モジュールを作成
以下のコードをコピペ
VBAのウインドウを消す
Alt +F8でVBAを起動
実行

色の変更シリーズ

選択セルの文頭・文末の色を変更

文字数の調整：c.Characters([文頭から何文字目], [文字数])

セルの文頭から

Sub 選択セルの文頭から色を変更()
For Each c In Selection
c.Characters(1, 1).Font.Color = RGB(255, 0, 0)
Next c
End Sub

セルの文末から

Sub 選択セルの文末から色を変更()
For Each c In Selection
c.Characters(Len(c), 1).Font.Color = RGB(255, 0, 0)
Next c
End Sub

選択セル中で検索して色を変更

こちらのページで紹介されている方法が便利です

vbabeginner.net

「Call 指定文字列のフォント変更("hoge", RGB(255, 0, 0), False)」のFalseをTrueにするとボールド体
複数の文字を一括変更：「Call 指定文字列のフォント変更("hoge", RGB(255, 0, 0), False)」

以下、自分用に改変。

Sub 指定文字列のフォント変更(a_sSearch, a_lColor, a_bBold)
    Dim f   As Font     '// Fontオブジェクト
    Dim i               '// 引数文字列のセルの位置
    Dim iLen            '// 引数文字列の文字数
    Dim r   As Range    '// セル範囲の１セル
    
    iLen = Len(a_sSearch)
    i = 1
    
    '// 選択セル範囲を１セルずつループ
    For Each r In Selection
        '// 指定されたセルの文字列から引数文字列を全て検索
        Do
            '// セル文字列から引数文字列を検索
            i = InStr(i, r.Value, a_sSearch)
            
            '// 引数文字列が存在しない場合
            If (i = 0) Then
                '// 次検索用に検索開始位置を1に初期化
                i = 1
                
                '// このセルの処理を終了
                Exit Do
            End If
            
            '// 引数文字列部分のFontオブジェクトを取得
            Set f = r.Characters(i, iLen).Font
            
            '// フォント設定
            f.Color = a_lColor  '// 文字色
            f.Bold = a_bBold    '// 太さ
            
            '// 次検索用に検索開始位置をずらす
            i = i + 1
        Loop
    Next
End Sub

Sub 指定文字列のフォント変更実行()
    Call 指定文字列のフォント変更("hoge", RGB(255, 0, 0), False)
End Sub

選択セル中で色文字の色を一括変更

www.excel.studio-kazu.jp

こちらのページに記載のものをいつも使用させていただいています。

Font.ColorIndex＝3　が赤
Font.ColorIndex =5　が青

 Sub 色変換()
    Dim i As Integer
    Dim c As Range
    For Each c In Selection    '選択された範囲内で
        For i = 1 To Len(c)
            With c.Characters(i, 1)
                If .Font.ColorIndex = 3 Then 'その文字が「赤」のときは
                    .Font.ColorIndex = 5     'その文字を「青」にする
                End If
            End With
        Next i
    Next c
 End Sub

追記シリーズ

置換で何とかなることのほうが多い。

選択セルに文字を追加

セルの文末に

Sub 選択したセルの文末に文字列追加()
　Selection.Value = Evaluate(Selection.Address & "&"" hogehoge""")
End Sub

セルの文頭に

Sub 選択したセルの文頭に文字列追加()
 Selection.Value = Evaluate("""hogehoge""&" & Selection.Address)
End Sub

2021-03-01

機械学習の概要

データサイエンス

機械学習の概要

機械学習の定義

機械学習の定義は以下の2つが有名です。

1つ目は、Arthur Samuel（ア－サー・サミュエル）による定義で次のようなものです。

The field of study that gives computers the ability to learn without being explicitly programmed.
（訳：明示的にプログラムしなくても、学習する能力をコンピュータに与える研究分野）

2つ目は、Tom Mitchell（トム・ミッチェル）の著書『Machine Learning』(1997) の序文にあるもので、次のような見解です。

A computer program is said to learn from experience E with respect to some class of tasks T and performance measure P, if its performance at tasks in T, as measured by P, improves with experience E.
（訳：コンピュータプログラムは､タスクTとパフォーマンス測定Pに関連する経験Eから学習すると言われています。タスクTでパフォーマンスをした場合、パフォーマンス測定Pによって達成度が測定され、経験Eによって改善されていきます。）

タスクTとは、「どのように課題を実行するか」であり、学習済みモデルが与えるソリューションを指します。その内容（＝機械学習でできること）は、回帰、分類、主成分分析、クラスタリング、異常検知、機械翻訳、転写、ノイズ除去などがあげられます。あくまで、モデルが期待される機能に関することであって、学習中に得られることとは切り離された概念です。例えば、「アヤメの画像から種類を分類する」というモデルを作る際には、分類自体がタスクＴなのであって、エッジ検出などは含まれません。

性能指標Pとは、「性能をどう評価するか」を指します。これにより、学習を行う指針を得ることができます（例：モデルが分類した結果が、どのくらい正確であったか）。

経験Eとは、データ・特徴量との関連を得られているかを指します。正解ラベルの有無（教師あり学習／教師なし学習）、画像、上の例のアヤメであれば、さらに茎の長さ、葉の大きさ、発芽時期、開花時期などが考えられるでしょう。

パラメトリックモデルとノンパラメトリックモデル

機械学習のモデルは、背景にある分布などを仮定するかどうかにより、次の二つに分類することができます。

パラメトリックモデル：あらかじめデータの背景にある確率分布を仮定し、パラメータだけを学習します。例えば、物理学や化学などの知見から「予測値は〇〇関数に従う」ということや、統計学の知見から「予測値は〇〇分布に従う」ということが予測できた場合、その関数に対してフィッティングを行うことなどがあります。このような前提を置くことにより、少ないデータで精度よく計算することができるようになります。一方、ある特定の条件だけに使用できる機械学習のモデルになってしまうデメリットもあります。

ノンパラメトリックモデル：データの背景にある分布や関数そのものを推定する方法です。物理学や化学、統計学の深い洞察が無くとも使用することができ、多くの事例に機械学習のモデルを適用できる可能性があります（高い汎用性）。しかし、パラメトリックモデルよりも多くのデータが必要となるデメリットもあります。運用面では機械の予測が「ブラックボックス化」してしまうと、せっかく構築したモデルが採用されないこともあるかもしれません。

というように分けられます。

教師あり学習と教師なし学習

教師あり学習（Surpervised Learning）とは、人間が正解データを作り、これを訓練データとして機械に教え込むこむ方法です。機械がある値を予測する際に、そのターゲットとなる値を目的変数と言います。それが他のデータ（説明変数）にどのように依存しているかを探ることが機械の役割になります。

教師なし学習（Unsupervised Learning）とは、人間が正解データを作るのではなく、機械自身がデータの背景にあるパターンや示唆を見つけてゆく手法です。データには「目的変数」はありません。機械は、データをいくつかのグループに分類したり（クラスタリング）、データの情報を失わないようにより小さい次元に縮約する主成分分析（Principle Component Analysis, PCA）などの手法があります。

強化学習

強化学習（Reinforcement Learning）とは、機械の行動（の戦略や結果）に対して報酬を与えることで、機械にベストプラクティスを見つけさせようとする手法です。強化学習は、「行動」に対する「報酬」（≒正解）があるという点では、教師あり学習に似ているといえます。逆に教師あり学習と異なる点は、部分ではなく全体で最適化を行うということです。例えば将棋で言えば、強化学習は「次の行一手」ではなく、数手先、あるいは勝負がつくまで計算を行い最適な行動を選択します。

マルチタスク学習

マルチタスク学習とは、単一のモデルで複数の課題を解く手法を指します。例えば、画像認識では領域の判定（セグメンテーション）とクラス分類、自然言語処理では品詞の分類、分節の判定、係り受けなどを同時に行うという事です。

バッチ学習とオンライン学習

機械学習のシステムはデータからモデル内にあるパラメータを適切に更新してゆくタイミングによっても分類できます。

バッチ学習では、全てのデータを利用して学習します。扱うデータ量にもよりますがそれなりに計算コストがかかるため、計算資源と時間（とお金）が必要になります。新しいデータがある場合は、ゼロから学習しなおし古いモデルと差し替えることになります。

オンライン学習データが追加されるたびに、パラメータの更新を行います。1データの学習のみなので、計算コストがかからず、その場ですぐにデータの学習をすますことができます。株価のように断続的にデータを受け取り、素早く変化に対応しなけれならないシステムに有効です。オンライン学習の弱点はデータにノイズなどがのってしまった場合、影響が大きく性能が下がってしまう事です。

ミニバッチ学習は、バッチ学習とオンライン学習の合いの子で、ある程度まとまったデータごとにパラメータの更新を行います。変化にそれなりに対応でき、ノイズにそれなりに強く、計算コストと時間（とお金）をそれなりに抑えることができます。

学習と検証

独立同分布仮定

機械学習のモデルには、未知のデータに対しても適切な結果を与えられることが求められます（汎化性能）。逆に、学習時に使用したデータでの（性能指標Ｐによる）性能は高いけれど、それ以外での性能が低い状態のことを過学習といいます。

そのため、通常、収集したデータを「訓練データ」と「テストデータ」に分け、「訓練データ」を用いて学習を行ったあとに「テストデータ」を用いて性能を評価するということを行います（ホールドアウト法）。

f:id:Shint:20210131212353p:plain

経験Ｅとなるデータに対して、 - 各データが互いに独立 - 訓練データとテストデータが同一の確率分布に従うという性質を独立同分布（Independent and Identically Distributed; IID）仮定と言います。

さらに、交差検証の中でもよく利用されるk分割交差検証法（Cross-Validation, CV）では、データをk個に分割して、そのうち1つをテストデータ、それ以外を訓練データとして性能を評価します。

f:id:Shint:20210131223749p:plain

過学習

学習を進めてゆくと「訓練データでは性能が高かったが、それ以外のデータではむしろ性能が低かった」ということが起こります。これを過学習といいます。訓練データに特化してチューニングを進めてしまったがために、未知のデータに関してはむしろ誤った答えを推測するという現象です。

下図のようなモデルでは、与えられた点に関しては精度よく再現していますが、少しでもパラメータ $\, x_i \,$ の値が異なると、途端にモデルの精度が下がるはずです。

f:id:Shint:20210214164542p:plain

正則化

過学習は、訓練データの量やノイズに対してモデルの方が複雑すぎるときに起こります。そこで、モデルに制約を加え、モデルを単純化することを正則化（regularization）といいます。

具体的なアルゴリズム

線形回帰

単純な回帰のモデルとして線形回帰があります。予測する値 $\, y \,$ が、入力するn個の特徴量 $\, \{ x_1, x_2, \cdots, x_n \} \,$ の一次関数で与えられると考えます。

$f(x_1, x_2, \cdots, x_n) = w_1 x_1 + w_2 x_2 + \cdots + w_n x_n + b \\$

とします。 $\, \{ w_1, w_2, \cdots, w_n \} \,$ は特徴量の係数で重みといい、また、 $\, b \,$ は定数項でバイアスといいます。

ここで、 $\, \boldsymbol{w} = (w_1, w_2, \cdots, w_n) \,$ および $\, \boldsymbol{x} = (x_1, x_2, \cdots, x_n) \,$ と定義して、

$f(\boldsymbol{x}) = w_1 x_1 + w_2 x_2 + \cdots + w_n x_n +b = \boldsymbol{w} \cdot \boldsymbol{x} + b\\$

のように表すと便利です。

例
あるコンビニのアイスクリームの販売個数 $\, y \,$ を予測するタスクを考える。アイスクリームは暑い日になれば売れると考えれば「気温」や「湿度」が特徴量になると考えられる。これらの変数を $\, x_1, x_2 \,$ と表記する $\, (n=2) \,$ 。

次に、販売個数 $\, y \,$ は $\, x_1, x_2 \,$ の一次関数であると仮定する。その際、「気温が1℃上昇」することと「湿度が1%上昇」することの売り上げに対する影響度は異なる。それゆえ、その影響度を個別に調整する重み $\, w_1, w_2 \,$ を導入して、

$y = w_1 x_1 + w_2 x_2 + b \\$

とする。ただし、 $\, b \,$ は定数項である。

さて、モデルが学習するというのは、この重み $\, \{ w_1, w_2, \cdots, w_n \} \,$ が手元にあるデータを再現するように調整するということです。予測値があるデータに近いかどうかを表す指標として平均二乗誤差（Root Means Square Error：RMSE）が一般的です。

$\displaystyle \begin{eqnarray*} RMSE &=& \sqrt{\dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i) } \\ &=& \sqrt{\dfrac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] } \\ \end{eqnarray*}$

$\, \{ y_1, y_2, \cdots, y_m \} \,$ は実際のデータで、 $\, m \,$ はデータの個数です。要は、この差を最小にするように重み $\{ w_1, w_2, \cdots, w_n \}$ を調整すればよいということになります。

なお、RMSEは標準偏差と似ていますが、RMSEは各データと推定値との差です（次のサイトが参考になります：二乗平均平方根などの統計量について）

例
この3日間の気温、湿度、アイスクリームの売り上げデータがあり、これらから明日のアイスクリームの販売個数を予想するタスクでは、

データ数： $m=3$
実際の販売個数： $y_1, \, y_2, \, y_3$
予測販売個数： $\displaystyle \, \boldsymbol{w} \cdot \boldsymbol{x}_1 +b, \, \boldsymbol{w} \cdot \boldsymbol{x}_2 + b, \, \boldsymbol{w} \cdot \boldsymbol{x}_3 +b \,$

RMSEを最小化する重み $\, \{ w_1, w_2, \cdots, w_n \} \,$ を考える際には、平方根を中身だけを考えればよいので、平均二乗誤差（Means Square Error：MSE）を扱うことが機械学習では一般的です。このように、学習を進めるにあたり最小化したい関数をコスト関数といいます。

$\displaystyle \begin{eqnarray*} MSE &=& \dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i) \\ &=& \dfrac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] \\ \end{eqnarray*}$

リッジ回帰

過学習は、訓練データの量やノイズに対してモデルの方が複雑すぎるときに起こり、モデルに制約を加えることで過学習を抑制するのでした（正則化）。

リッジ回帰（Ridge Regression）は線形回帰の正則化で、コスト関数に $\, \alpha/2 \sum_{i=1}^{n} w_i^{2} \,$ という正則化項を加えます。リッジ回帰のコスト関数を $\, J(\boldsymbol{w}) \,$ とすると、

$\displaystyle \begin{eqnarray*} J(\boldsymbol{w}) &=& \dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i) + \dfrac{\alpha}{2} \sum_{i=1}^{n} w_i^{2}\\ &=& \dfrac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] + \dfrac{\alpha}{2} \sum_{i=1}^{n} w_i^{2} \\ \end{eqnarray*}$

このようにすることで、学習モデルは訓練データを再現するだけでなく、重み $\, \{ w_1, w_2, \cdots, w_n \} \,$ できる限り小さく保つという束縛条件が一つ加わることになります。

ここで、 $\, \alpha \,$ は人間が設定する値（ハイパーパラメータ）です。 $\, \alpha =0 \,$ とすればモデルはただの線形回帰となり、 $\, \alpha \,$ を大きな値にすれば、重み $\, \{ w_1, w_2, \cdots, w_n \} \,$ はすべてがほとんどゼロになり、モデルの予測値はほぼ定数となります。

数学的には、正則化項としてベクトルの大きさを表すノルムを加えていることになります。リッジ回帰の場合は $\, L^{2} \,$ ノルムとなります。

参考
$\, n \,$ 次元のベクトル $\, \boldsymbol{x} = (x_1, x_2, \cdots, x_n) \,$ および、 $\, p ( 1\leq p \leq \infty ) \,$ に対して、 $\, \sqrt[p]{|x_1|^p+|x_2|^p+\cdots+|x_n|^p}$ を $\, \boldsymbol{x} \,$ の $\, L^p \,$ ノルムといい、 $\,|| \boldsymbol{x}||_p \,$ と表す。

ラッソ回帰

ラッソ回帰（Least Absolute Shrinkage and Selection Operator Regressio: Lasso Regression）はリッジ回帰と同様に、コスト関数に正則化項を加えますが、 $\, L^{2} \,$ ノルムではなく、 $\, L^{1} \,$ ノルムを使います。

$\displaystyle \begin{eqnarray*} J(\boldsymbol{w}) &=& \frac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i) + \alpha \sum_{i=1}^{n} |w_i| \\ &=& \frac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] + \alpha \sum_{i=1}^{n} |w_i| \\ \end{eqnarray*}$

ラッソ回帰は重要性の低い特徴量の重みを0にする傾向があります。言い換えれば、自動的に必要な特徴量を選択し疎なモデルとなるということです。

しかし、ラッソ回帰は特徴量の数が多い時や特徴量の間に強い相関があるときには挙動が不規則になることがあることが注意点です。

Elastic Net

ラッソ回帰の弱点を克服しようとするのがErasitc Netです。考え方は、リッジ回帰とラッソ回帰の合いの子です。

$\displaystyle \begin{eqnarray*} J(\boldsymbol{w}) &=& \dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i) + r \alpha{2} \sum_{i=1}^{n} |w_i| + \dfrac{1-r}{2} \alpha \sum_{i=1}^{n} |w_i|^{2} \\ \end{eqnarray*}$

$\, r=0 \,$ でリッジ回帰、 $\, r=1 \,$ でラッソ回帰と等しくなります。

リッジ回帰で試した際に意味がある特徴量が疑われる際には、ラッソ回帰の代わりにElastic Netを使用する方が良いとされています。

ロジスティック回帰

ロジスティック回帰（Logistic Regression）は、回帰という名前がついていますが、分類の用途で使われます。

確率 $\, p \,$ であるカテゴリに属すことの起こりやすさを表す量として、オッズ比 $\, p/(1-p) \,$ があります。そして、この対数をとった関数をロジット関数 $\, L(p) \,$ といいます。

$L(p) := \log \dfrac{p}{1-p} \\$

$\, L(p) \,$ は定義域 $\, 0 \leq p \leq 1 \,$ で値域 $\, -\infty \,$ ～ $\, \infty \,$ となる関数です。ここで、ロジット関数の値は、特徴量 $\, \{ x_1, x_2, \cdots, x_n \} \,$ の線形結合で与えられると仮定します。

まず、特徴量の線形関係を（線形回帰の時と同様に）

$\begin{eqnarray*} z &=& w_1 x_1 + w_2 x_2 + \cdots + w_n x_n + b \\ &=& \boldsymbol{w} \cdot \boldsymbol{x} + b \\ \end{eqnarray*}$

としますロジット関数がこのような特徴量の1次関数 $\, z \,$ で記載されているとすると、

$\log \dfrac{p}{1-p} = z \\$

これを $\, p \,$ について解くと、

$p = \frac{1}{1+e^{-z}} \\$

と、シグモイド関数となります。

実装としては、特徴量 $\, \{ x_1, x_2, \cdots, x_n \} \,$ の線形結合をシグモイド関数に代入すればよいわけですから、以下のようになります。

import numpy as np

def input(x, w, b):
    z = np.dot(x,w) + b
    y = 1/(1+np.exp(-z))
    return y

ロジスティック回帰の学習とは線形回帰同様に、重み $\, \boldsymbol{w} \,$ 、バイアス $\, b \,$ を適切に調整してゆくという事になります。それでは、コスト関数はどう考えるかというと、ここで確率統計で学んだ最尤法を用いてゆきます（確率統計を参照してください）。

シグモイド関数の形からわかるように、 $\, p \geq 0.5 \,$ ならば「あるカテゴリに属す」、 $\, p \lt 0.5 \,$ ならば「あるカテゴリに属さない」という分類ができます。つまり、二値分類となりますので、

$t_i =\begin{cases} 1 & (p \geq 0.5) \\ 0 & (p \lt 0.5) \\ \end{cases}$

となる $\, t_i (i=1,2,\cdots, n) \,$ を導入すると、i番目の $\, t_i \,$ が $\, 1 \,$ か $\, 0 \,$ になる確率は、 $\, p^{t_i} (1-p)^{1-t_i} \,$ とまとめてかけるので、n個分類した時に実現する同時確率 $\, L(\boldsymbol{w},b) \,$ は

$\displaystyle L(\boldsymbol{w},b) = \prod_{i=1}^{n} p_i^{t_i} (1-p_i)^{1-t_i} \\$

となります。この対数をとったものが対数尤度となります。ここでは、（最大問題を最小問題として扱うため）マイナスをかけた値をとり、

$\displaystyle \log L( \boldsymbol{w}, b ) = - \sum_{i=1}^{n} \{ t_i \log p_i + (1-t_i) \log (1-p_i) \} \\$

とします。この形の関数を交差エントロピー誤差（Cross Entropy Loss）関数と言います。

# クロスエントロピー誤差
def closs_entropy(p, t):
    loss = -np.sum( (t * np.log(p) + (1 - t) * np.log(1 - p)), axis=0)
    return loss

この関数を最小にするような $\, p \,$ が「もっともらしい」値となります。

ソフトマックス回帰（多項ロジスティック回帰）

ロジスティック回帰は2つのカテゴリへの分類でした。 3つ以上の分類ができるように、シグモイド関数をソフトマックス関数に差し替えたモデルを、 ソフトマックス回帰といいます。

$\displaystyle p_i := \dfrac{e^{z_i}}{\sum_{j=1}^{n} e^{z_j} } \quad where \quad z_i = \boldsymbol{w_i} \cdot \boldsymbol{x} +b_i \\$

# ソフトマックス回帰
def softmax(x, w, b):
    z = np.dot(w, x) + b
    p = np.exp(z) / np.sum(np.exp(z))
    return p

k近傍法

k近傍法（k-Nearest Neighbors: k-NN）とは、既知のグループが複数ある時に、どのグループに属するかわからない新たなデータが、どのグループに属するのかを推定する手法です。

f:id:Shint:20210215230723p:plain

上図はk近傍法のイメージです。赤と青で示したデータがあった際に、新しいサンプルがどちらのグループに属するかを考えます。データの空間において新しいサンプルから近いデータを3つ選ぶと $\, (k=3) \,$ 、多数決により新しいサンプルは青のグループに分類されます。

k近傍法のハイパーパラメタは次の2つです。

距離
- 数学的に「距離の公理」を満たす定義は多数あります（ユークリッド距離、マハラノビス距離など）。
新しいサンプルの近傍点の個数
- 上の図で $\, k=3 \,$ では青に分類されますが、 $\, k=7 \,$ では赤に分類されるように結果が異なり得ます。

上の例は分類問題の例でしたが、回帰にも使えます。

分類問題の場合：近傍k個のサンプルのうち最も多いクラスを予測クラスとする
回帰問題の場合：予測値として近傍k個のサンプルの平均値をとる

k近傍法の特徴としては、

ノンパラメトリックモデル
オンライン学習が容易：パラメータの更新がないため、新しいデータを加えるだけで予測が可能です

サポートベクトルマシン

サポートベクトルマシン（Support Vector Machine）の考え方を以下に示します。

f:id:Shint:20210222154524p:plain

二つのクラスに対し境界線を与える際に、サポートベクトルマシンは境界線と隣接する点の距離（マージン）が最大となるように境界線を決定します。この隣接するデータ点のことをサポートベクトルと言います（あるデータ点 $\, (x_0, x_1, \cdots, x_n) \,$ は一つのn次元ベクトルとしてみることができるため「ベクトル」です）。

線形分離可能なデータに関して：ハードマージン

下図のように線形分離可能な場合の定式化をハードマージンといい、以下、2次元の例を用いてマージンはどのように最大化されるか詳しく見てゆきます（n次元でも大筋は同様です）。

赤のサポートベクトルで示される点 $\, S_{+} \,$ の位置を表すベクトルを $\, \boldsymbol{X}_{S+} = ( X_{1S+}, X_{2S+})^{t} \,$ 、そこから境界線 $\, l \,$ におろした垂線の足 $\, H_{+} \,$ の位置を表すベクトル $\, \boldsymbol{X}_{H+} = (X_{1H+}, X_{2H+})^{t} \,$ とします。

まず、境界線 $\, l \,$ を

$\begin{eqnarray*} w_1x_{1}+w_2x_{2} + b =0 \qquad \cdots \quad (1) \\ \end{eqnarray*}$

とすると、 $\, H_{+} \,$ に対しては、

$\begin{eqnarray*} w_1X_{1H+}+w_2X_{2H+} + b =0 \qquad \cdots \quad (2) \\ \end{eqnarray*}$

という関係が成り立ちます。

次に、境界線 $\, l \,$ に対する直行するベクトルを考えます。オレンジには、境界線 $\, l \,$ の係数を成分とするベクトル $\, \boldsymbol{w} \,$ を示しました。一方、境界線 $\, l \,$ は、 $x_{1} \,$ が $\, w_{2} \,$ だけ進むと、 $x_{2} \,$ が $\, w_{1} \,$ だけ減少する関係を示したものですから（図中の水色）、図を見比べると、 $\, \boldsymbol{w} \,$ は境界線 $\, l \,$ と直行するベクトルということがわかります（大学の数学的には、これは(1)の左辺のgradientをとることでも容易に得られます）。

今、 $\, \vec{H_+S_+} \,$ と $\, \boldsymbol{w} \,$ は平行であるため、

$\begin{pmatrix} X_{1S+} - X_{1H+} \\ X_{2S+} - X_{2H+} \end{pmatrix} = t \begin{pmatrix} w_1\\ w_2 \end{pmatrix} \qquad \cdots \quad (3) \\$

ただし、tはゼロでない定数です。以下、 $\, \vec{H_+S_+} \,$ の大きさを求めるために、 $\, t \,$ を計算します。

式(3)の両辺で $\, \boldsymbol{w} \,$ との内積をとり

$w_1(X_{1S+}-X_{1H+})+w_2(X_{2S+}-X_{2H+}) =t(w_1^2+w_2^2) \\$

展開して、

$\begin{eqnarray*} w_1X_{1S+}+w_2X_{1S+} + b \\ - (w_1X_{1H+}+w_2X_{2H+} + b) \end{eqnarray*} =t(w_1^2+w_2^2) \\$

左辺のカッコ内は式 $\, (2) \,$ より0となるため、

$\ t = \dfrac{w_1X_{1S+}+w_2X_{2S+} + b}{w_1^2+w_2^2} = \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{S+} + b }{|\boldsymbol{w}|^2} \\$

よって、

$\vec{H_+S_+} = \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{S+} + b }{|\boldsymbol{w}|} \begin{pmatrix} w_1/|\boldsymbol{w}| \\ w_2/|\boldsymbol{w}| \end{pmatrix} \\$

$\, \vec{H_+S_+} \,$ の大きさがサポートベクトルマシンにおけるマージンで、それを $\, M \,$ とします。また、Class赤に分類されるサポートベクトル以外のデータ点と境界線 $\, l \,$ の距離は $\, M \,$ より大きくなければいけませんから、 $\, x_{i'} \,$ をClass赤に分類されるデータ点として、

$\dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i'} + b }{|\boldsymbol{w}|} \geq M \\$

が成立します。等号はサポートベクトルの点において成立します。この関係性を保ちながら $\, M \,$ を最大化することが、今、得たい結果となります。境界式 $\, l \,$ を動かす際のパラメタは、式(1)で表されるように $\, \{ w_1, w_2, b \} \,$ または $\, \{ \boldsymbol{w}, b \} \,$ となりますので、

$\displaystyle \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i'} + b }{|\boldsymbol{w}|} \geq M, \quad \max_{\boldsymbol{w},b } M \qquad \cdots \quad (4) \\$

f:id:Shint:20210223170847p:plain

一方、クラス青に分類されるサポートベクトル $\, S_- \,$ から境界線 $\, l \,$ におろした垂線の足 $\, H_- (X_{1H-}, X_{2H-}) \,$ とで得られるベクトル $\, \vec{H_{-}S_{-}} \,$ は、境界線 $\, l \,$ と反平行になりますので、符号が逆になります。

$\displaystyle - \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i''} + b }{|\boldsymbol{w}|} \geq M, \quad \max_{\boldsymbol{w},b } M \qquad \cdots \quad (5) \\$

ここに、 $\, x_{i''} \,$ はClass青に分類されるデータ点です。

Class赤では正、Class青では負の符号を表す $\, y_i \,$ を導入すると、式(4)(5)はまとめて次のように書くことができます。

$\displaystyle y_i \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i} + b }{|\boldsymbol{w}|} \geq M, \quad \max_{\boldsymbol{w},b } M \qquad \cdots \quad (6) \\$

ここで、 $\, x_i \,$ は任意のデータ点です。

境界線 $\, l \,$ の両辺を定数倍しても結果は変わらないため、計算しやすいように、以下のような変換を行います。

$\displaystyle \tilde {\boldsymbol{w}} = \dfrac{\boldsymbol{w}}{M | w | }, \quad \tilde {b} = \dfrac{b}{M | w | } \\$

この変換の式。および、変換後のベクトルの大きさ $\, | \tilde{w} | = 1/M \,$ より、

$\displaystyle y_i ( \tilde{ \boldsymbol{w} } \cdot \boldsymbol{X}_{i} + \tilde{b} ) \geq 1, \quad \max_{ \tilde{ \boldsymbol{ w }}, b } \dfrac{1}{ | \tilde{w} | } \\$

2つめの条件について、最大値をとることは逆数の最小値をとることと等価ですから、

$\displaystyle y_i ( \tilde{ \boldsymbol{w} } \cdot \boldsymbol{X}_{i} + \tilde{b} ) \geq 1, \quad \min_{ \tilde{ \boldsymbol{ w }}, b } \dfrac{ | \tilde{w} |^2 }{2} \\$

というようにして、計算をします（2乗にしたり、1/2を掛けたりするのは、最適化計算をしやすくしているだけなので、気にしないでOKです）。

これらの関数を最適化してゆくアルゴリズムに受け渡すことで学習が進みます。

線形分離不可能なデータに関して：ソフトマージン

これまでの例は線形分離可能な場合を扱ってきましたが、そうでない場合（線形分離不可能な場合）の定式化をソフトマージンといいます。

f:id:Shint:20210301211219p:plain

図のように、境界線から片側のマージン分だけ離れたところにできる境界を超平面といいます（ここでの例は「直線」ですが、3次元に一般化すると「平面」、3次元以上では「超平面」と呼びます）。今は、Class 赤側を正（Class 青側を負）としていました。

線形分離不可能な場合、各超平面を超えたデータ（分離できないデータ）に対してペナルティ（損失）を与えます。

損失の与え方は幾通りも考えられますが、一番自然なのは上の図のように「超平面を超えたデータに対して、超平面からの距離を損失としてあたえる」かと思います。

具体的に考えると、

あるデータ $\, X_i \,$ が超平面を超えない場合（ $\, y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b) \leq M \,$ ）、ペナルティは0
あるデータ $\, X_i \,$ が超平面を超える場合（ $\, y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b) \gt M \,$ ）、ペナルティは超平面からの距離

f:id:Shint:20210301210907p:plain — 図．Class 赤に対するペナルティの与え方の例

となるので、まとめて、

$\displaystyle \xi_i = \max \left\{ 0, \, M - \dfrac{y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b)}{||w||} \, \right\} \\$

のように書けます。結局、ソフトマージンSVMの行うことは、「マージン $\, M \,$ を最大にしつつペナルティ合計を最小化するような境界線 $\, l \,$ を与える」ということになります。

「最大」と「最小」が混在していますので、少々めんどくさいですね。これは、ハードマージンの時と同じように変数変換をすれば、最小化問題に集約できます。つまり、変数変換後のペナルティ $\, \tilde{\xi}_i \,$

$\displaystyle \tilde{\xi}_i = \max \left\{ 0, \, 1 - y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b) \, \right\} \\$

として、

$\displaystyle y_i ( \tilde{ \boldsymbol{w} } \cdot \boldsymbol{X}_{i} + \tilde{b} ) \geq 1 - \tilde{ \xi}_i, \quad \min_{ \tilde{ \boldsymbol{ w }}, b } \left\{ \dfrac{ | \tilde{w} |^2 }{2} + C \sum_{i=1}^{n} \tilde{\xi}_i \right\} \\$

となります。ここであらたに登場した $\, C (\leq 0)\,$ は、ハイパーパラメタで二つの量のバランスを調整します。

以上の議論に出てきたペナルティ項をスラック変数といいます。

2021-02-12

文字の被りを回避してワード文章の一括置換する方法

データサイエンス

背景

目的

日本語で記載した設計書をもとに、日本語変数名とコーディングに使う変数名（VarName）を追記して、検索の手間を省きたい！と考えました。

今回は、エクセルで「”日本語変数名", "日本語変数名【VarName】"」のような置換テーブルを作成し、このテーブルをもとに一括置換するプログラムを作成します。

課題

今回課題となったのは、長い日本語変数名のなかに他の短い変数名が含まれてしまっていたことでした。

f:id:Shint:20210212215846p:plain — 図．長い変数名に短い変数名が含まれるイメージ

上の「式」単位で配列に保存されているイメージです。

今回作ったもの

内容

今回は、簡易的に「長い変数名を選択する」というロジックを組んでみました。

f:id:Shint:20210212221136p:plain — 図．日本語変数名が長い候補で置換するイメージ

一つの配列要素のなかに、短い変数名とそれを含む長い変数名が無かったことは救いでした。次のようなものがあると、短い方が取り残されます。

f:id:Shint:20210212221648p:plain — 図．今回のロジックでは同一の式に短い変数名があるとNG

あと、最大の日本語変数名の長さが2つある場合、および、同一の日本語変数名を使いまわしている場合は、その変数だけ変換されません。

ソース

import docx
import time
import pandas as pd

print('The strat time is:', time.ctime())


def GetReplaceCSVtoList(path, filename):
    data = pd.read_csv(path + filename + '.csv', encoding='shift-jis')
    data = data.values.tolist()
    return data


def GetWordData(path, filename):
    doc = docx.Document(path+filename+'.docx')
    return doc


def ReplaceMaxLength(line, data):
    # cell内でのreplace候補を抽出
    temp = []
    for i in range(len(data)):
        if data[i][0] in line:
            temp.append(data[i])

    # cell内のreplace候補のうち、他の変数名の一部である変数を抽出
    temp2 = []
    for i in range(len(temp)):
        temp1 = temp[0:i]+temp[i+1:len(temp)]  # 自分以外のリスト
        temp1 = [word[0] for word in temp1]  # 自分以外のリストのkey
        for k in range(len(temp1)):
            if temp[i][0] in temp1[k]:
                temp2.append(temp[i][0])

    # cell内でのreplace候補から、他の変数名の一部である変数を削除
    temp3 = [word for word in temp if not (word[0] in temp2)]

    # 変換
    for i in range(len(temp3)):
        line = line.replace(temp3[i][0], temp3[i][1])

    return line


# --- main ---
path = './（変換テーブルを記載したCSVを含んだフォルダ名）/'
filename = '（ファイル名）'
data = GetReplaceCSVtoList(path, filename)
data = [[word[2], word[3]] for word in data]

path = './（処理前のワードが保存されているフォルダ名）/'
filename = '（処理前のワード名）'
doc = GetWordData(path, filename)

flag = '（処理する表の一つ目のヘッダー）'
for tbl in doc.tables:
    rowtemp = tbl.rows[0]
    if rowtemp.cells[0].text == flag:
        for row in tbl.rows:
            for cell in row.cells:
                cell.text = ReplaceMaxLength(cell.text, data)

doc.save('（保存するワード名）.docx')

print('The finish time is:', time.ctime())

2021-02-12

wordファイルから特定の表の内容を抽出する方法

データサイエンス

目的

背景

前回、二つの.txtファイルの各列を比較し、相手側の.txtファイルに含まれていない文字列を返してくれるスクリプトを作成しました。

shint.hatenablog.com

目的

今回の目的は、「wordファイルに埋め込まれた表の中の文章を抽出し、比較する」です。

作ったもの

「ワードに埋め込まれた表から文章を抽出」の部分をUPします。手元では、前回のモノとマージして使っています。

動作環境

Python 3.7.4
python-docx 0.8.10

コード

#モジュール名はpython-docx(インストール時に注意)
import docx

doc = docx.Document("./（wordファイル名）.docx")
print("段落の個数：", len(doc.paragraphs))
print("表の個数：", len(doc.tables))
print("図（行内）の個数：", len(doc.inline_shapes))

allLine = []
for i in range(len(doc.tables)):
    tbl = doc.tables[i]
    row = tbl.rows[0]
    if row.cells[0].text == '（抽出したい表のヘッダーを記入してください）':
        for row in tbl.rows:
            values = [cell.text for cell in row.cells]
            allLine.append(values)
    else:
        continue

allLine

前回のものとマージしたもの

前回はリストを比較していましたが(数式の両辺で表記ゆれをチェック - Shint’s blog)今回はword内の表（2列）の比較です。自分用のモノなので、ソースにメモ書きしちゃっています。

# ワードファイルが保管されているフォルダ名を記入してください。
word_path = './WordData/'

#「WordData」（上記で記入したフォルダ名）フォルダに対象のワードファイルを置いてください。
# Wordの機能で「校閲/承諾/すべての変更を反映し、変更の記録を停止」を行って保存してください。

# ワードファイル名を記入してください（例『hoge.docx』）
word_filename = 'hoge'

# チェックするワードファイル中の表のヘッダー1列目を記入してください。
# これ以外のヘッダーはskipされます
flag = '日本語変数名'

# まず、「＋+×-−－–÷／＝=Σ{}｛｝><＞＜≧≦・\n」でリストが分割されます。
# その上で不要な要素を『Stopword.txt』に記載してください（1行、1要素）
# その上で要素内の不要な語句を『Stemmer.txt』に記載してください（1行、語句）
# 『Stopword.txt』、『Stemmer.txt』などが含まれるフォルダ名を記入してください
inp_path = './inp/'


# 結果が出力される「res」フォルダを準備してください。

# 左辺のまとめが「LHS.txt」、右辺のまとめが「RHS.txt」に出力されます。
# 左辺で定義されていない右辺の日本語変数名が「Results.txt」に出力されます。

######################################################

import docx
import re
import itertools
import time

# ファイルを開く
def OpenFile(path, filename):
    file = open(path+filename+'.txt', 'r', encoding='utf-8')
    data = file.read()
    data = data.split('\n')
    return data


# 式の前後で分割
def Splitter(templist):
    for i in range(len(templist)):
        templist[i] = re.split('[＋+×-−－–÷／＝=Σ{}｛｝><＞＜≧≦・\n]', templist[i])
    templist = list(itertools.chain.from_iterable(templist))
    return templist


# StopWord
def StopWord(templist):
    templist = [templist[i].strip()
                for i in range(len(templist))
                if not(templist[i].strip() in StopWordList or templist[i].strip() == '')
                ]
    return templist


# 不要な文字列を削除(Stemmer)
def Stemmer(templist):
    templist1 = []
    for i in range(len(templist)):
        for k in range(len(StemmerList)):
            templist[i] = templist[i].replace(StemmerList[k], '')
        templist1.append(templist[i].strip())
    templist1 = [line.strip() for line in templist1 if not(line.strip() == '')]
    templist1 = list(dict.fromkeys(templist1))
    return templist1


# 2つのリストの比較
def compare(lisname1, listname2):
    templist = [lisname1[i]
                for i in range(len(lisname1))
                if not((lisname1[i] in listname2))
                ]
    return templist


# ワードから直接リストを取得
def GetTableData(path, filename, flag):
    doc = docx.Document(path+filename+'.docx')
    print('【', filename+'.docxの情報 】')
    print('    段落の個数：', len(doc.paragraphs))
    print('    表の個数：', len(doc.tables))
    print('    図（行内）の個数：', len(doc.inline_shapes))
    allline = []
    for i in range(len(doc.tables)):
        tbl = doc.tables[i]
        row = tbl.rows[0]
        if row.cells[0].text == flag:
            for row in tbl.rows:
                values = [cell.text for cell in row.cells]
                allline.append(values)
        else:
            continue
    return allline


# 2×2のリストから列を取得
def GetColumn(data, ColumnNum):
    temp = []
    for i in range(len(data)):
        if len(data[i]) != 2:
            print('Error: 表は2列で記載してください')
        elif data[i][0] == data[i][1]:
            continue
        else:
            temp.append(data[i][ColumnNum])
    return temp


# 前処理
def ProcessingColumn(templist):
    templist = Splitter(templist)
    templist = StopWord(templist)
    templist = Stemmer(templist)
    return templist

# ----------------------------------------------------

print('Calculation started at', time.ctime())

# ワードファイルから表を抽出
data = GetTableData(word_path, word_filename, flag)


# 表から列を抽出
column0 = GetColumn(data, 0)
column1 = GetColumn(data, 1)
data = []

# 入力データを読み込み
parameters = OpenFile(inp_path, 'parameter')
StopWordList = OpenFile(inp_path, 'StopWord')
StemmerList = OpenFile(inp_path, 'Stemmer')

# 右辺前処理前を保存
l_out = open('./res/RHS_before.txt', 'w', encoding='utf-8')
for i in range(0, len(column1)):
    print(column1[i], file=l_out)

# 前処理
column0 = ProcessingColumn(column0)
column1 = ProcessingColumn(column1)
parameters = ProcessingColumn(parameters)

# 右辺前処理後を保存
f_out = open('./res/RHS.txt', 'w', encoding='utf-8')
for i in range(0, len(column1)):
    print(column1[i], file=f_out)

# 左辺前処理後を保存
g_out = open('./res/LHS.txt', 'w', encoding='utf-8')
for i in range(0, len(column0)):
    print(column0[i], file=g_out)

# 2つのリストを比較し、未定義変数を出力
result = compare(column1, parameters) #パラメータリスト分を除去
result = compare(result, column0)

h_out = open('./res/Results.txt', 'w', encoding='utf-8')
for i in range(0, len(result)):
    print(result[i], file=h_out)

f_out.close()
g_out.close()
h_out.close()

print('Calculation finished at', time.ctime())

2021-02-08

数式の両辺で表記ゆれをチェック

データサイエンス

目的

背景

設計書などを作る際に、WordやExcelなどで、

新変数名 ＝パラメタ × 既存の変数名

と記載することがある。

当然ながら、既存変数名はそれまでに左辺で定義していなければならない。さらに、それらは同一の名前でなければならない（パラメタは別ファイルで入力値として与えるものとする）。

式の本数が少ないアプリケーションの場合は問題にならないが、すべてマニュアルで記載する場合、1000本の数式を超えると記載ミス（半角/全角、スペース有/無などといった変数名の揺らぎ）も発生する。

やりたいこと・課題

右辺のパラメタと既存変数名が左辺の追加変数名で定義されているかチェックしたい。

その際の課題としては、以下のようなものがある。

wordで変更の履歴を記録していると検索に引っかかりにくい（目視確認時には「シンプルな変更履歴/コメント」にすることを推奨）
wordでは半角が改行位置にあると隠れてしまい判別が不能。
一方、wordではあいまい検索ができてしまうため、その機能を停止しないで進めるとゆらぎが気づきにくい（マニュアルで検索する場合は、「高度な検索」から「あいまい検索」のチェックボックスを外すことを推奨）。
「Σ」や「～の場合」など、簡単なチェックツールを使用するには不要な文字がある。

今回作ったもの

内容

左辺と右辺のゆらぎのチェックツールを作成した。

四則演算の記号などでリストを分割
不要なリストの要素を削除をできるようにした（Stop Word）
不要なリストの要素の一部を削除をできるようにした（Stemmer）
左辺と右辺をそれぞれ.txtファイルで与えて、上の処理を行った後に比較し、定義されていない文字を返す。
パラメタも.txtファイルで与えることを想定。

使い方

フォルダ「データ」にデータファイルを保存。
数式の左辺を「colmun1.txt」に改行区切りで保存（この際、不要な行、文字があっても可）。
数式の右辺を「colmun2.txt」に改行区切りで保存（この際、不要な行、文字があっても可）。
パラメタを「parameters.txt」に改行区切りで保存。
数式にif文などでコメントがあることを想定。例えば『(比較したい式) < 0』から『(比較したい式)』を抽出するため、数学記号で分割。何で分割するかはSplitter()内に記入する。
不要な配列要素を「StopWord.txt」に記入。
- 計算を走らせて、不要なところをコピペする。
配列要素の一部分の不要な語句を「Stemmer.txt」に記入。
- 計算を走らせて、不要なところをコピペする。
比較したいリストを本文のcompare(lisname1, listname2)に記入。
- listname1からlistname2を引いたリストが「res」フォルダに出力される。

動作環境

Python 3.7.4

コード

import re
import itertools


# ファイルを開く
def OpenFile(path, filename):
    file = open(path+filename+'.txt', 'r', encoding='utf-8')
    data = file.read()
    data = data.split('\n')
    return data


#式の前後で分割
def Splitter(templist):
    for i in range(len(templist)):
        templist[i] = re.split('[＋×−÷=Σ><≧≦]', templist[i])
    templist = list(itertools.chain.from_iterable(templist))
    return templist


#StopWord
def StopWord(templist):
    templist = [templist[i].strip()
                for i in range(len(templist))
                if not(templist[i].strip() in StopWordList or templist[i].strip() == '')
                ]
    return templist


# 不要な文字列を削除(Stemmer)
def Stemmer(templist):
    templist1 = []
    for i in range(len(templist)):
        for k in range(len(StemmerList)):
            templist[i] = templist[i].replace(StemmerList[k], '')
        templist1.append(templist[i].strip())
    templist1 = [line.strip() for line in templist1 if not(line.strip() == '')]
    templist1 = list(dict.fromkeys(templist1))
    return templist1


# 前処理
def ProcessingColumn(templist):
    templist = Splitter(templist)
    templist = StopWord(templist)
    templist = Stemmer(templist)
    return templist

# 2つのリストの比較
def compare(lisname1, listname2):
    templist = [lisname1[i]
                for i in range(len(lisname1))
                if not((lisname1[i] in listname2))
                ]
    return templist

#--- main ---

# データの読み込み
path = './データ/'
column1 = OpenFile(path, 'column1')
column2 = OpenFile(path, 'column2')
StopWordList = OpenFile(path, 'StopWord')
StemmerList = OpenFile(path, 'Stemmer')
parameters = OpenFile(path, 'parameter')


# 前処理
column1 = ProcessingColumn(column1)
column2 = ProcessingColumn(column2)


#比較
result = compare(column2, parameters) #パラメタ分を除去
result = compare(result, column1)


#出力
f_out = open('./res/results.txt', 'w', encoding='utf-8')
for i in range(0, len(result)):
    print(result[i], file=f_out)

2021-02-08

Pythonチートシート（データ加工）

データサイエンス

これは自分用のチートシートで随時更新型エントリです

Import

とりあえず書いておく

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from jupyterthemes import jtplot
%matplotlib inline
jtplot.style(theme='solarizedd')
import pickle

たまに使う

import seaborn as sns
import re
import os
from tqdm._tqdm_notebook
import tqdm_notebook
tqdm_notebook.pandas(desc='progress')
import time

ロガー

import logging
logging.basicConfig(filename='log.txt', level=logging.DEBUG, format='%(levelname)s:%(message)s')
ogger = logging.getLogger(__name__)

データ処理の基本

CSVからの読み込み

data = pd.read_csv('../00_data/name.csv', encoding='SHIFT-JIS')
display(data)

ファイル内の名前を一括変更

import glob
import os

path='C:\\Users\\huga\\hoge\\'
filename = glob.glob('C:\\Users\\huga\\hoge\\*')

for i in range(0, len(filename)):
    filename0 = filename[i].split('\\')[-1]
    filename1 = filename[i].split('\\')[-1].replace('_English', '')
    os.rename(path+filename0, path+filename1)

縦方向に結合(concat)

concat_data = pd.concat([data1,data2], ignore_index=True)

横方向に結合(join)

huga = pd.merge(data1[['id', 'column1', 'column2']],
         data2[['id', 'column1', 'column2']],
         on = 'id',
         how = 'left')

data1とdata2でkeyの書き方が異なる場合（例：「customer」と「顧客名」など）

huga = pd.merge(data1[['id', 'column1', 'column2']],
         data2[['id', 'column1', 'column2']],
         left_on = 'id_left',
         right_on = 'id_right',
         how = 'left')

グループにまとめる

join_data.groupby(['まとめる列名']).count()

join_data.groupby(['row','column']).sum()[['value1', 'value2']]

pivotテーブル

join_data.groupby(['row','column']).sum()[['value1', 'value2']]
pd.pivot_table(join_data,
               index='row',
               columns='column',
               values=['value1','value2'],
               aggfunc='sum')

df.pivot_table(index='row',
               columns='column',
               aggfunc='size',
               fill_value=0)

複数のリストの中から、合計値が大きいリストを選択

#例えばa1~a3のようなリストがあるとする
a1 = [1,2,3,4,5]
a2 = [6,7,8,9,10]
a3 = [11,12,13,14,15]

a = [a1, a2, a3]
a_sum = [sum(x) for x in a]
a = [np.argmax(a_sum)]

空の要素を削除

templist = [line.strip() for line in templist1 if line.strip() != '']

データ記入の揺れ補正

#ユニークな数データ数がいくつあるか確認（あるべき数より多い＝データに揺らぎがある）
print(len(pd.unique(df['column_name'])))
#大文字・小文字、全角・半角スペースの補正
df['column_name'] = df['column_name'].str.upper() #全て大文字に
df['column_name'] = df['column_name'].str.replace('　','') #全角スペースを削除
df['column_name'] = df['column_name'].str.replace(' ','')  #半角スペースを削除
#データ名でソートすることによりデータを見やすくする
df.sort_values(by=['column_name'], ascending=True)
#ユニークなデーを確認
print(pd.unique(df['column_name']))
#ユニークなデータの数があるべき数であるかを確認
print(len(pd.unique(df['column_name'])))

欠損値を他から補完

flg_is_null = df['column_name1'].isnull() 
for trg in list(df.loc[flg_is_null, 'column_name2'].unique()):
    #nullがあるdf['column_name1']を選び、nullでないところからdf['column_name2']を抽出する    
    value_column_name1 = df.loc[
        (~flg_is_null)&(df['column_name2']==trg), 
        'column_name2'
    ].max()

    #nullがある商品名を選び、nullでないところから商品価格を抽出する
    df['column_name1'].loc[
        (flg_is_null)&(df['column_name2']==trg)
    ] = value_column_name1
display(df.head())
display(df.isnull().sum()) #すべての値がゼロであることを確認

loc[<インデックスラベル>, <カラムラベル>] というようにインデックスラベルにTtue/Falseを入れるとTrueの部分だけだけ取り出される
.unique()をつけることで、'column_name1'が欠損している'column_name2'の重複をなく抽出する。
'~'は否定演算子で、flg_is_null==Falseと同じ
&は比較演算で、比較するものが同じなら1、違うのならば0。
1==TrueはTrue
今回はmaxで一つの値に絞り込み。本来は、「2つの値があったらエラー」等のロジックがあってもよい。
出した値が正しいかは次のロジックでもよい

for trg in list(df['column_name2'].sort_values().unique()):
    print(
        trg
        + 'の最大額：' 
        + str(df.loc[df['column_name2']==trg]['column_name1'].max())
        + 'の最小値：'
        + str(df.loc[df['column_name2']==trg]['column_name1'].min(skipna=False))
    )

日付の形式の統一

数値データを日付に変換

flg_is_serial = df['date'].astype('str').str.isdigit()
fromSerial = pd.to_timedelta(data.loc[flg_is_serial, 'data'].\
    astype(float), unit="D") + pd.to_datetime('1900/01/01')
fromSerial

.str.isdigit()は「全ての文字が数字なら真、そうでなければ偽」 to_timedeltaでDatetime型に変換。

次に、念のため日付データも変換

fromString = pd.to_datetime(data.loc[~flg_is_serial, 'data'])
fromString

変換したデータを結合

df['date'] = pd.concat([fromSerial, fromString])

保存・ログ

Pickle

#保存
with open('picklename.binaryfile', 'wb') as web:
    pickle.dump(to_pickle, web)

#ロード
with open('picklename.binaryfile', 'rb') as web:
    from_pickle = pickle.load(web)

ロガー

logger.debug(f'(名前): {(引数)}')

演算子

代数演算子

通常の「+,-,*,/」以外に

a % b         # a を b で割った余り
a ** b        # a の b 乗
a // b        # 切り捨て除算

ビット演算子

~a            # ビット反転(1と0を反転)
a & b         # AND:論理積(aもbも1のビットが1)
a | b         # OR:論理和(aまたはbが1のビットが1)
a ^ b         # XOR:排他的論理和(aまたはbが1のビットが1)
a << b        # b ビット左シフト
a >> b        # b ビット右シフト

2021-01-19

確率統計

データサイエンス

有名なディープラーニングの資格にJDLAのE資格があります。

本エントリーは必要な確率統計の知識のミニマムを解説します。

確率統計
確率分布の例
最尤法
参考書

確率統計

期待値・分散・標準偏差

期待値、分散、標準偏差の定義式を以下に示します。

離散確率変数Xが確率関数 $\, P_i (i=1,2, ･･･, n) \,$ の確率分布に従うとき、
　(1) 期待値 $\mu = E [ X ] = \sum_{i=1}^{n} x_i P_i$
　(2) 分散 $V [ X ] = \sum_{i=1}^{n} (x_i - \mu) P_i = E [ (X- \mu )^{2} ]$
　(3) 標準偏差 $\sigma = \sqrt{V [ X ] }$

(2)については、以下の式が成立します。

　 $\sigma^{2} =V [ X ] = E [ X^{2} ] - E [ X ]^2$

証明

$\begin{eqnarray*} V [ X ] &=& E [ (X - \mu)^{2} ] \\ &=& E [ X^{2} -2 \mu X + \mu^{2} ] \\ &=& E [ X^{2} ] - 2 \mu E [ X ] +E [ X ]^{2} \\ &=& E [ X^{2} ] - E [ X ]^2 \end{eqnarray*}$

これらは連続変数の場合にも成り立ちます。

連続確率変数Xが確率密度 $\, p(x) \,$ の確率分布に従うとき、
　(1) 期待値 $\, \mu = E [ X ] = \int_{\infty}^{-\infty} x p(x) dx$
　(2) 分散 $\, V [ X ] = \int_{\infty}^{-\infty} (x - \mu)^{2} p(x) dx = E [ X^{2} ] - E [ X ]^{2}$
　(3) 標準偏差 $\, \sigma = \sqrt{V [ X ] }$

ネイマン・ピアソン統計

いわゆる普通の統計学のことをネイマン・ピアソン統計と言います。データの持つ性質などから適切な統計分布を用いて表現し、その分布に関する命題について真偽の判断を与えようとする手続きです（統計的仮説検証）。

統計分布の表現方法として2つの種類があります。1つ目は記述統計学。これはデータを整理して、表やグラフなどを用いてそのデータの傾向や特徴をとらえることを言います。

2つ目は推測統計学です。全体の中から抽出されたサンプル（標本）から全体（母集団）の傾向や特徴を推測します。全体の平均がどのくらいであるかなど推測を行ったり、2つのデータ間の平均値に有意差があるかなどといった仮説を検定したりします。

具体的な手法をゼロから知りたい方には、『完全独習統計学入門』（小島寛之著）をお勧めします（が、E資格という観点では不要かと思います）。

ベイズ統計

ベイズ統計（とネイマン・ピアソン統計とベイズ統計の違い）を分かりやすく説明している本として、『完全独習　ベイズ統計学入門』（小島寛之著）がすこぶる楽しいです。

ベイズ統計は、新しい情報が得られたときにその前後でどのように確率が変化するかを与えます。それは、条件付確率

$P(A|B)$

という形で表されます。Aは事象、Bは前提条件を表します（つまり、「ある前提条件Bのもとで、ある事象Aが起こる確率」という意味です）。

例えば、 $\, x, y \,$ という2つの事象について（小島先生に従って）次の図を用いて説明します。(a)は2つの事象が同時に起こる確率を表しています。左上の領域が「 $\, x \,$ と $\, y \,$ が同時に起こる確率」でオレンジの枠が全体です。よって面積比としては、 $\, x\,$ と $\, y \,$ が同時に起こる確率 $\, P(x,y)=1/4 \,$ となります。

仮に、ここで「 $\, y\,$ が起こった」という情報があったとします（前提条件）。すると、「 $\, y \,$ が起こらなかった」領域は考えなくてよいので、面積比として条件付き確率 $\, (P(x|y)=)1/2 \,$ というように変化します。

f:id:Shint:20210123225632p:plain

例題
東京の人口を1000万人として、新型コロナウィルスに罹患している人（感染力を持つ人）が10万人いたと仮定する。今回、東京の全人口にPCR検査を行うこととなり、あなたがPCR陽性になったとする。この場合、本当にあなたが罹患している確率はどの程度だろうか？ただし、PCR検査の精度は以下のものとする。

コロナウイルスに感染している人に対して、70%の割合で陽性反応を示す。
コロナウイルスに感染していない人に対しも、1%の割合で陽性反応を示す。

解説

感染者／非感染者の人数、陽性／陰性（偽陽性／偽陰性）を図にまとめると次のようになります。

f:id:Shint:20210124224434p:plain

ここに、 (a) 10万人×0.7、(b)10万人×0.3、(c)990万人×0.01、(d)990万人×0.99 です。今、あたなたは、「PCR陽性である」ということが確定したわけですから、図の黄色でハッチングしたところのみ考えればよいわけです。

つまり、東京で全数検査をした際に自分が陽性になってしまったときに、本当に「自分は感染者かどうか」という確率は、

$\frac{10^5 \times 0.7}{10^5 \times 0.7+9.9\times 10^6 \times 0.01} =41.4(\%)$

となります。意外に低い値ではないでしょうか？全数検査の費用対効果を考えると、仮に感染力のある人数が10万人いても、なかなか実施するのは難しく感じます。真の陽性者 (a) 10万人×0.7 = 7万人に対し、自宅待機ならともかく、真面目に対策をする場合、偽陽性者(c)990万人×0.01=9.9万の隔離病棟を準備しなくてはいけません。

ベイズの定理

条件付き確率 $\, P(y|x) \,$ がわかっているときに $\, P(x|y) \,$ が知りたいときに便利なのがベイズの定理です。

$P(x|y) = \frac{P(x)P(y|x)}{P(y)} \\$

分母は、 $\, P(y) = \sum_{x} P(y|x)P(x) \,$ で計算できます。

例題
転職を望むあるグループの中で、JDLAの資格を持った人が以下のように分類できるとする。

事象A：E資格（およびG資格）を保有
事象B：G資格のみ保有
事象C：E資格もG資格も保有していない

さらに、〇×テクノロジー株式会社への転職活動の結果を次のように分類する。

事象E：〇×テクノロジー株式会社の中途採用に合格
事象F：〇×テクノロジー株式会社の中途採用に不合格

これらの割合は、以下のように与えられるとする。

$P(A) = 0.1 \qquad P(B) = 0.5 \qquad P(C) = 0.4 \\$

$P(E|A) = 0.9 \qquad P(E|B) = 0.6 \qquad P(E|C) = 0.2 \\$

この時、〇×テクノロジーに合格できた人が、E資格に合格できる割合を求めよ。

解説

ベイズの定理より、

$\begin{eqnarray*} P(A|E) &=& \frac{P(A)P(E|A)}{P(A)P(E|A)+P(B)P(E|B)+P(B)P(E|B)} \\ &=& \frac{0.1 \times 0.9}{0.1 \times 0.9 + 0.5 \times 0.6 + 0.4 \times 0.2} \\ &=& 0.191 \end{eqnarray*}$

というように求めることができます。先ほどの図のように表すと、

f:id:Shint:20210124233614p:plain

となり、(a)0.1×0.9、(b)0.5×0.6、(c)0.4×0.2、(d)0.1×0.1、(e)0.5×0.4、(f)0.4×0.6です。

ここで、「〇×テクノロジーに合格した」という追加の情報があれば、不合格の部分は考えなくてよくなるのですから、黄色のハッチング部分のみを考えればよいことになります。黄色のハッチング部分を全体として、(a)の領域の割合を出すことが $\, P(A|E) \,$ を求めるということになります。

確率分布の例

ベルヌーイ分布

何かが起こった時に、起こる結果が「A or B」のように2つに分類できる試行のことベルヌーイ試行と言います。例えば、サイコロを振った時に「事象A：5以上の目」「事象B：4以下の目」のようなケースです。

Aとなる確率を $\, p \,$ 、Bとなる確率を $\, q \,$ とし、

$P(A) = p \qquad P(B) = q = 1 - p$

と表すことにします。

このようなベルヌーイ試行をn回行い、ちょうどk回だけ事象Aとなる確率は次のように表されます。

f:id:Shint:20210125221514p:plain

①は「 $\, n \,$ 回中 $\, k \,$ 回選択する組み合わせ」の意味でコンビネーション、②は「確率 $\, p \,$ が $\, k \,$ 回」、③は「確率 $\, q \,$ が $\, n-k \,$ 回」を意味します。このように、確率変数 $X=k$ に対して( $\, q=1-p \,$ を用いて)

$P(X) = {}_n \mathrm{C} _k q^k (1-p)^{n-k} \quad (k=0,1,2, \cdots, n)$

であらわされる分布をベルヌーイ分布または二項分布といいます。

ベルヌーイ分布では次の関係性が成り立ちます。
　(1) 期待値 $E [ X ] = np$
　(2) 分散 $V [ X ] = npq$

(1)の証明

コンビネーションを階乗であらわし、 $\, n \,$ と $\, p \,$ で括り、二項定理を使用します。また $\, p+q=1 \,$ の関係も使います。

$\begin{eqnarray*} E [ X ] &=& \sum_{k=0}^n k \cdot {}_n \mathrm{C}_k \cdot p^k q^{n-k} \\ &=& \sum_{k=0}^n k \frac{n!}{k!(n-k)!} \times p^k q^{n-k} \\ &=& np \sum_{k=1}^n \frac{(n-1)!}{(k-1)!(n-k)!} \times p^{k-1} q^{n-k} \\ &=& np \sum_{k=0}^{n-1} \frac{(n-1)!}{k!(n-k-1)!} \times p^{k-1} q^{n-k-1} \\ &=& np \sum_{k=0}^{n-1} {}_{n-1} \mathrm{C}_k \cdot p^{k-1} q^{n-k-1} \\ &=& np \cdot (p+q)^{n-1} \\ &=& np \end{eqnarray*}$

(2)の証明

分散 $\, V [ X ] \,$ の関係式より、

$\begin{eqnarray*} V [ X ] &=& E [ X^{2} ] - E [ X ]^{2} \\ &=& E [ X^{2} ] - E [ X ] + E [ X ] - E [ X ]^{2} \\ &=& E [ X^{2} - X ] + E [ X ] - E [ X ]^{2} \\ &=& E [ X(X-1) ] + E [ X ] - E [ X ]^{2} \qquad \cdots \quad (1) \\ \end{eqnarray*}$

となりますから、まず、右辺第1項 $\, E [ X(X-1) ] \,$ を計算します。上の式と同じ方針で、 $\, n(n-1)p^{2} \,$ で括り、残りの部分には二項定理を適用します。

$\begin{eqnarray*} E [ X(X-1) ] &=& \sum_{k=0}^n k(k-1) \cdot {}_n \mathrm{C}_k \cdot p^k q^{n-k} \\ &=& \sum_{k=0}^n k(k-1) \frac{n!}{k!(n-k)!} \times p^k q^{n-k} \\ &=& n(n-1)p^2 \sum_{k=2}^n \frac{(n-2)!}{(k-2)!(n-k)!} \times p^{k-2} q^{n-k} \\ &=& n(n-1)p^2 \sum_{k=0}^{n-2} \frac{(n-2)!}{k!(n-k-2)!} \times p^{k-2} q^{n-k-2} \\ &=& n(n-1)p^2 \sum_{k=0}^{n-2} {}_{n-2} \mathrm{C}_k \cdot p^{k-2} q^{n-k-2} \\ &=& n(n-1)p^2 \cdot (p+q)^{n-2} \\ &=& n(n-1)p^2 \\ \end{eqnarray*}$

この結果を(1)に代入して、

$\begin{eqnarray*} V [ X ] &=& n(n-1)p^2 + np - (np)^2 \\ &=& np(np-p+1-np) \\ &=& np(1-p) \end{eqnarray*}$

E資格では $\, n=1 \,$ のケースが問われるようです。

ベルヌーイ分布
　確率変数 $X( x = 0,1)$ について、
　　(1) 確率分布 $P(X) = p^x (1-p)^{1-x}$
　　(2) 期待値 $E [ X ] = p$
　　(3) 分散 $V [ X ] = p(1-p)$

マルチヌーイ分布

ベルヌーイ分布が「事象A or 事象B」のように2つに分類できる問題に対して利用できる分布であったのに対し、 3つ以上に分類できる問題にも対応できるように拡張した分布がマルチヌーイ分布です。

いずれか1つのみ「 $\, 1 \,$ 」の値をとり、それ以外「 $\, 0 \,$ 」の値をとる $\, n \,$ 個のパラメタ{ $\, x_1, x_2, \cdots, x_n \,$ }と、それに対応する確率{ $\, p_1, p_2, \cdots, p_n \,$ } ( $\, \sum_{k=1}^n p_k = 1 \,$ )について

$P(X=x) = \prod_{k=1}^N p_k^{x_k}$

となる確率分布です。

例えば、サイコロの目「3」のみが $\, X=1 \,$ それ以外が $\, 0 \,$ 、それぞれの目が出る確率は $\, 1/6 \,$ とすれば、 Xは、マルチヌーイ分布に従います。

正規分布

最もよく使われる分布に正規分布（またはガウス分布）があります。

$N(x; \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2\sigma^2} \right)$

$\, \mu \,$ は平均、 $\, \sigma^{2} \,$ は分散( $\, > 0 \,$ )、 $\, \sigma \,$ は標準偏差を表します。

$\, \mu = 0 \,$ 、 $\, \sigma^{2} = 1 \,$ である場合、標準正規分布といいます。

$N(x; 0, 1) = \frac{1}{\sqrt{2 \pi }} \exp \left( -\frac{x^2}{2} \right)$

任意の正規分布は標準正規分布に変換することができます。確率変数 $\, X \,$ が正規分布 $\, N(\mu, \sigma^{2}) \,$ に従うとき、上の2つの式（の特に指数）を見比べて、変数変換 $\, Z=\frac{X-\mu}{\sigma^{2}} \,$ をすれば、 $\, N(0, 1) \,$ となることがわかるかと思います。

最尤法

「現実となった確率変数は、確率が最大のものが実現した」という仮定を最尤原理（principle of maximum likelihood）といいます。

例えば、まず、「１」となる確率が $\, p \,$ 、「０」となる確率が $\, 1-p \,$ のベルヌーイ分布があり、 $\, ( x_1, x_2, x_3, x_4, x_5 ) = (1,1,1,1,0) \,$ の $\, n=5 \,$ の確率変数があるとしましょう。

確率Lは、

$\displaystyle L(p) = p^4(1-p)$

であらわされるわけですから、試しに $\, p = 0.1 \,$ と $\, p = 0.9 \,$ の時を考えてみますと、

$L(p=0.2) = 0.1^{4} \times 0.9 = 9 \times 10^{-5} \,$
$\, L(p=0.8) = 0.9^{4} \times 0.1 = 0.0651\,$

となり、後者の方が確率 $\, L(p=0.9) \,$ が大きく「もっともらしい」といえます（今の場合、計算するまでもなく、4回出てくる「1」の確率の方が大きいはずです）。

一般化すると、 $\, L(p) \,$ は確率 $\, p \,$ のとりうる値におけるもっともらしさを表す関数とみなすことができ、このような「もっともらしさ」を尤度（likelihood）、その関数を尤度関数（likelihood function）といいます。

最尤法（maximum likelihood method）とは、尤度関数を最大にする $\, p \,$ を考えることとなります。

上の例では、

$\displaystyle L(p) = p^3(4-5p) \\$

なので $\, p=0.8 \,$ が最大値となります。5回中4回が「1」となるので、 $\, p=4/5 \,$ と推測できることに一致しています。

一般に、確率分布の変数を $\, \theta \,$ として、n個の確率変数 $\, \{ x_1, x_2, \cdots, x_n \} \,$ の同時確率分布は $\, \theta \,$ の関数とみなし $\, f(x, \theta) \,$ とすると、尤度関数は確率関数の積となり、

$\displaystyle L(\theta) = f(x_1, \theta) \, f(x_2, \theta) \, \cdots \, f(x_n,\theta) = \prod_{i=1}^{n} f(x_n, \theta) \\$

となります。ただし、正規分布が（ $\, \mu, \sigma^{2} \,$ ）の関数であるように、 $\, \theta \,$ は必ずしも1つではありません。

尤度関数の形として、積よりも対数をとって和の形にした方が扱いやすいため、

$\displaystyle \log L(\theta) = \sum_{i=1}^{n} \log f(x_i, \theta) \\$

を考えるのが普通です（対数尤度）。現実に起こりそうな $\, \theta \,$ とは、対数尤度を最大にする $\, \theta \,$ ですから、

$\displaystyle \frac{\partial \log L(\theta)}{\partial \theta} = 0 \\$

Yahooショッピング

7net