MeCabの辞書の生起コストの調整で分割誤りを直す

以前、MeCabのユーザー辞書に単語を追加する方法を紹介しました。
参考: MeCabでユーザー辞書を作って単語を追加する
(この記事の本題とずれますが、M1/M2チップのMacではいろんなファイルのパスが変わっているのでこちらの記事を参照される際はご注意ください。)

実際、先ほどの記事を参照してどんどん単語を追加してくと、追加したのに形態素解析結果に反映されないということがあります。というよりも、そもそもデフォルトのIPA辞書でさえ正しい単語が追加されているのにそれが使われずに誤った分割がされることがあります。

このような場合、多くのケースで単語の生起コストを修正することで対応可能です。
シンプルに、MeCabに使って欲しい単語の生起コストを下げるだけです。ただ、あまりに極端に下げてしまうと逆にその低すぎる生起コストのせいで分割すべき場合でも分割されないという事象が発生しえます。

そこで、どの程度生起コストを下げたらいい感じで分割ができるようになるのかというのを探っていこうというのがこの記事です。

例として、今回は「高等学校」って単語でテストしていきましょう。いろんな個別の高等学校の学校名が固有名詞としてIPA辞書に登録されているのに、なぜか一般名詞の「高等学校」が登録されてないのですよね。

ユーザー辞書の登録の記事に沿ってseedファイルを用意し、コストを推定すると、
高等学校,1285,1285,5078,名詞,一般,,,,,高等学校,コウトウガッコウ,コートーガッコー
となり、生起コストは5078となります。

これをコンパイルして、 sample.dic ってユーザー辞書を作っておきます。

さて、適当に例文を作って形態素解析してみましょう。

import pandas as pd
import subprocess
import MeCab


# サンプルデータ生成
sample_texts = [
    "九州高等学校ゴルフ選手権",
    "地元の高等学校に進学した",
    "帝京高等学校のエースとして活躍",
    "開成高等学校117人が現役合格",
    "マンガを高等学校の授業で使う",
]
df = pd.DataFrame({"text": sample_texts})

# MeCabの辞書のパス取得
dicdir = subprocess.run(
    ["mecab-config", "--dicdir"],
    capture_output=True,
    text=True,
).stdout.strip()
sysdic = f"{dicdir}/ipadic"

# 分かち書き用のTagger
tagger_w = MeCab.Tagger(f"-O wakati -d {sysdic} -u sample.dic")

print(df["text"].apply(tagger_w.parse).str.strip())
"""
0          九州 高等 学校 ゴルフ 選手権
1       地元 の 高等 学校 に 進学 し た
2     帝京 高等 学校 の エース として 活躍
3    開成 高等 学校 117 人 が 現役 合格
4     マンガ を 高等 学校 の 授業 で 使う
Name: text, dtype: object
"""

はい、以上の最後の結果の通り、「高校」と「学校」の間にスペースが入っていてせっかく追加した「高等学校」は使われませんでした。ちなみに、Nベスト解を出すと2番目以降に登場するので、ユーザー辞書には正常に追加されています。

print(df["text"].apply(lambda x: tagger_w.parseNBest(2, x)).str.strip())
"""
0                九州 高等 学校 ゴルフ 選手権 \n九州 高等学校 ゴルフ 選手権
1          地元 の 高等 学校 に 進学 し た \n地元 の 高等学校 に 進学 し た
2      帝京 高等 学校 の エース として 活躍 \n帝京 高等学校 の エース として 活躍
3    開成 高等 学校 117 人 が 現役 合格 \n開成 高等学校 117 人 が 現役 合格
4      マンガ を 高等 学校 の 授業 で 使う \nマンガ を 高等学校 の 授業 で 使う
"""

なぜこうなるのか、というと「高等」の生起コストが924、「学校」の生起コストが1624とどちらもとても低いんですよね。両単語間の連接コストが1028ありますが、それを考慮しても、924+1624+1028=3576で、「高等学校」一単語の生起コスト5078よりも低いです。

また、「高等」は名詞,形容動詞語幹なので、前の単語との連接コストも変わってきます。

そのため、追加した「高等学校」を単語として使ってもらうためには、生起コストをモデルが推定した値そのままではなくもっと低く設定してあげる必要があります。

さて、ここからどのくらい生起コストを下げていくらにしたらいいのかを探っていきましょう。方法としては、前後の品詞の単語との組み合わせをありうる全パターン考慮して生起コスト、連接コストから計算することも可能ですが、これは極端に低い値に設定する必要が出ることもあり、逆に分割誤りを誘発することもあるので、もっと具体的に正しく分割されて欲しい例文を基準に計算するのがおすすめです。この記事では先ほど使った5文を使います。

その方法は単純で、「高等」「学校」に割れてしまう分割のコストの総和と、制約付き解析で求める「高等学校」を強制的に使った分割のコストの差分を求めてその分を調整します。

# 生起コスト+連接コストの累積だけを返すTagger
tagger_pc = MeCab.Tagger(f"-F '' -E %pc -d {sysdic} -u sample.dic")
# それの制約付き解析版
tagger_p_pc = MeCab.Tagger(f"-p -F '' -E %pc -d {sysdic} -u sample.dic")

target_word = "高等学校"

# コストの合計をそれぞれ求める
df["コスト合計"] = df["text"].apply(lambda x: tagger_pc.parse(x)).astype(int)
# 制約付きの方は、制約をつけたい単語は 単語\t* 形式に置き換えて実行
df["制約付き解析コスト合計"] = df["text"].apply(
    lambda x: tagger_p_pc.parse(x.replace(target_word, f"\n{target_word}\t*\n"))
).astype(int)
# コストの差分を計算
df["コスト差分"] = df["制約付き解析コスト合計"] - df["コスト合計"]

# 結果表示
print(df)
"""
              text  コスト合計  制約付き解析コスト合計  コスト差分
0     九州高等学校ゴルフ選手権   6901         7495    594
1     地元の高等学校に進学した  10811        11807    996
2  帝京高等学校のエースとして活躍  15034        15484    450
3  開成高等学校117人が現役合格  40123        40138     15
4   マンガを高等学校の授業で使う  10407        12945   2538
"""

最後の例だけコスト差分がとにかく多いのは、を(助詞,格助詞,一般)と高等(名詞,形容動詞語幹)の連接コストが-1550と非常に小さいのも影響しています。

さて、これらの5文章で分割を正確にしたかったら、「高等学校」生起コストを2539下げれば大丈夫そうです。ただ、その前にちょっと実験として、996だけ小さくして、$5078-996=4082$に設定してみましょう。0,2,3番目はこれで正常に修正されるはず、4番目の文は修正されないはず、1番目はどうなるかみてみたいですね。

高等学校,1285,1285,4082,名詞,一般,,,,,高等学校,コウトウガッコウ,コートーガッコー

でユーザー辞書 sample2.dic を作ってみました。実験結果がこちら。

# 分かち書き用のTagger
tagger_w_2 = MeCab.Tagger(f"-O wakati -d {sysdic} -u sample2.dic")

print(df["text"].apply(tagger_w_2.parse).str.strip())
"""
0          九州 高等学校 ゴルフ 選手権
1      地元 の 高等 学校 に 進学 し た
2     帝京 高等学校 の エース として 活躍
3    開成 高等学校 117 人 が 現役 合格
4    マンガ を 高等 学校 の 授業 で 使う
Name: text, dtype: object
"""

ほぼ想定通りで、1番目のやつはシステム辞書の単語たちがそのまま使われましたね。もう1だけ下げる必要があったようです。

結果は省略しますが、もう一つだけコストを下げて4081にすると、1番目の文も正しく「高等学校」が単語として扱われます。

$5078-2538-1=2539$まで下げれば5文とも正しくなりますね。元のコストから見るとかなり低くしたように見えますが、「学校」の生起コストが1624なでのまぁ許容範囲でしょう。

以上の処理をメソッド化して、適切なコストを返してくれる関数を作るとしたらこんな感じでしょうか。あまり綺麗なコードじゃないのでもう少し考えたいですね。

wordで渡す引数はすでにユーザー辞書に登録済みで品詞等は判明していてあとはコスト調整だけの問題になっていることや、sentenceの中では正確に1回だけ登場することなど制約が多いのでかなり利用に注意が必要です。頻繁に必要になるようならまた改めて考えようと思います。

dicdir = subprocess.run(
    ["mecab-config", "--dicdir"],
    capture_output=True,
    text=True,
).stdout.strip()
sysdic = f"{dicdir}/ipadic"
# 単語の生起コストを出すTagger
tagger_p_c = MeCab.Tagger(f"-p -F %c -E '' -d {sysdic} -u sample.dic")
# コストの総和を出すTagger
tagger_pc = MeCab.Tagger(f"-F '' -E %pc -d {sysdic} -u sample.dic")
# 制約付き解析でコストの総和を出すTagger
tagger_p_pc = MeCab.Tagger(f"-p -F '' -E %pc -d {sysdic} -u sample.dic")


def calc_adjust_cost(word, sentence):
    all_cost = int(tagger_pc.parse(sentence))
    all_p_cost = int(tagger_p_pc.parse(sentence.replace(word, f"\n{word}\t*\n")))
    current_cost = int(tagger_p_c.parse(f"\n{word}\t*\n"))
    return current_cost - all_p_cost + all_cost - 1


calc_adjust_cost("高等学校", "マンガを高等学校の授業で使う")
# 2539

結果自体は正しそうですね。

PyGWalkerでデータフレームを可視化してみる

最近、TwitterでPyGWalkerというライブラリを知りました。
Githubで最初のコミットを見ると、2023/2/16の登場でまだ1ヶ月もたってない新しいライブラリです。

これが何かというと、JupyterでTableau風にDataFrameを可視化するライブラリになります。
(あくまでもTableau風であって互換というわけではないですが。)
Tableau好きな自分にとってはキャプチャ画像で見た第一印象が衝撃的だったので早速使ってみることにしました。

PyPiのページはこちらです。
参考: pygwalker · PyPI

インストールは普通にpipでできます。この記事を書いてる時点で入ったバージョンは pygwalker==0.1.4.9.post0 でした。

$ pip install pygwalker

早速これ使ってみましょう。使い方は簡単で、import して、walkメソッドに可視化したいデータフレームを渡すだけです。
とりあえず、scikit-learnのサンプルデータでやってみましょうかね。特に理由はないですが特徴量が比較的多様なwineを使います。まずデータフレーム作成から。

import pandas as pd
from sklearn.datasets import load_wine


# データ読み込み
wine = load_wine()
# 特徴量名を列名としてDataFrame作成
df = pd.DataFrame(
    wine.data,
    columns=wine.feature_names,
)

# target列も作成する。
df["target"] = wine.target
df["class"] = df["target"].apply(lambda x: wine["target_names"][x])

そして、次のコードでPyGWalkerが起動できます。Githubのサンプルに沿ってpygの名前でインポートしました。

import pygwalker as pyg


gwalker = pyg.walk(df)

そうすると以下のUIが表示されます。だいぶTableauそっくりじゃないですか?

とりあえず特徴量を2個選んで散布図でも作ってみましょう。prolineを X-Axis, flavanoidsをY-Axisのところにドラッグします。また、classで色分けしたいので、ColorのところにClassを持っていきます。この時点で、classごとに集計された3点の散布図になってるので、上の画像で黒反転してる立方体をクリックして集計を解除し個別プロットに変形します。グラフが少し小さいので、Layoutの設定もいじって拡大します。結果できるのが次のグラフです。

非集約化の捜査がちょっと独特でしたが、だいぶTableauライクな操作でいい感じの図が作れたと思います。

色を好みのものに変えたりとか、軸の細かい設定とかそういうTableauでならできる機能が色々なかったりしますが、それは今後のバージョンアップに期待しましょう。

オレンジのアイコンでグラフの種類も変えられます。折れ線グラフ、棒グラフ、円グラフなどの定番に加えて、Boxプロットとかもあるのもいいですね。
たとえば、Boxプロットを使うと次のようなグラフが作れます。

今のままの使い勝手だったら、to_csv して本家Tableauで表示するか、matplotlib使うかする方が便利かなぁとも思うのですが、なんせ登場してまだ1ヶ月未満で、バージョン0.1台のライブラリであり、頻繁に更新が入っている状況です。
これからの発展に期待したいですね。

Numpyでコレスキー分解と修正コレスキー分解を行う

最近書いたSVARモデルの話や、以前書いたVARモデルの直交化インパルス応答の話と少し関係あるのですが、コレスキー分解と言う行列分解があります。これをNumpyで行う方法を紹介しておこうと言う記事です。

参考: コレスキー分解 – Wikipedia

Wikipediaでは正定値エルミート行列について定義が書いてありますが一旦、実行列をメインで考えるのでその直後に書いてある実対称行列の話で進めます。
正定値実対称行列$\mathbf{A}$を下三角行列$\mathbf{L}$とその転置行列$\mathbf{L}^T$の積に分解するのがコレスキー分解です。

NumpyやScipyには専用のメソッドが実装されおり、容易に行うことができます。NumpyでできるならNumpyでやってしまうのがオススメです。(scipyの方はなぜか上三角行列がデフォルトですし。)
参考: numpy.linalg.cholesky — NumPy v1.24 Manual
scipy.linalg.cholesky — SciPy v1.10.1 Manual

では適当な行列用意してやってみましょう。

import numpy as np


# 適当に正定値対称行列を用意する
A = np.array([
    [2, 1, 3],
    [1, 4, 2],
    [3, 2, 6],
])

# 一応正定値なのを見ておく
print(np.linalg.eigvals(A))
# array([8.67447336, 0.39312789, 2.93239875])

# コレスキー分解とその結果を表示
L = np.linalg.cholesky(A)
print(L)
"""
[[1.41421356 0.         0.        ]
 [0.70710678 1.87082869 0.        ]
 [2.12132034 0.26726124 1.19522861]]
"""

# 転置行列と掛け合わせることで元の行列になることを確認
print(L@L.T)
"""
[[2. 1. 3.]
 [1. 4. 2.]
 [3. 2. 6.]]
"""

簡単でしたね。

さて、ここまでだとメソッド一つ紹介しただけなのですが、実はここからが今回の記事の本題です。Wikipediaには修正コレスキー分解ってのが紹介されています。これは元の行列を対角成分が全て1の下三角行列$\mathbf{P}$と、対角行列$\mathbf{D}$と、1個目の三角行列の転置行列$\mathbf{P^T}$の3つの積に分解する方法です。(Wikipedia中の式ではこの三角行列もLになってますが、先出のLと被ると以降の説明が書きにくいのでここではPとします。)
$$\mathbf{A}=\mathbf{PDP^T}$$
と分解します。

直交化インパルス応答関数の記事中で使ってるのはこちらなのですが、どうやらこれをやってくれる関数が用意されてないのでやってみました。
いやいや、このメソッドがあるじゃないか、と言うのをご存知の人がいらっしゃいましたら教えてください。

さて、僕がやった方法ですがこれは単純で、コレスキー分解の結果をさらに分解し、
$$\mathbf{L}=\mathbf{PD^{1/2}}$$
とすることを目指します。

ちょっと計算するとすぐわかるのですが、下三角行列$\mathbf{L}$を対角成分が1の下三角行列と対角行列の積に分解しようとしたら、対角行列の方は元の三角行列の対角成分を取り出したものにするしかありません。あとはそれ使って$\mathbf{P}$も計算したら完成です。

対角成分を取り出した行列はnumpy.diagを使います。これ2次元配列に使うと対角成分を1次元配列で取り出し、1次元配列を渡すとそこから対角行列を生成すると言うなかなかトリッキーな動きをするので、対角成分を取り出して対角行列を作りたかったらこれを2回作用させます。平方根も外すために2乗することも忘れないようにすると、$\mathbf{D}$は次のようにして求まり、さらにその逆行列を使って$\mathbf{P}$も出せます。

# 対角行列Dの計算
D = np.diag(np.diag(L))**2
print(D)
"""
[[2.         0.         0.        ]
 [0.         3.5        0.        ]
 [0.         0.         1.42857143]]
"""

# Dを用いて、対角成分1の下三角行列Pを求める
P = L@np.linalg.inv(D**0.5)
print(P)
"""
[[1.         0.         0.        ]
 [0.5        1.         0.        ]
 [1.5        0.14285714 1.        ]]
"""

# PDP^tで元の行列Aが復元できることを見る
print(P@D@P.T)
"""
[[2. 1. 3.]
 [1. 4. 2.]
 [3. 2. 6.]]
"""

はい、これでできました。ブログの都合上、逐一printして出力も表示してるのでめんどくさそうな感じになってますが、実質以下の3行で実装できています。

L = np.linalg.cholesky(A)
D = np.diag(np.diag(L))**2
P = L@np.linalg.inv(D**0.5)

statsmodelsでSVARモデルの推定

前回に続いてSVARモデル(構造VARモデル)の話です。前回はモデルの数式の形を紹介しただけだったので、今回は実際に仮想的なデータを作ってみてPythonのコードで推定してみます。

サンプルとして、以下のようなモデルを考えました。

$${\scriptsize
\left[\begin{matrix}1&0&0\\0.3&1&0\\-0.2&-0.3&1\end{matrix}\right]\mathbf{y_t}=
\left[\begin{matrix}0.8\\-0.5\\0.3\end{matrix}\right]+
\left[\begin{matrix}-0.1&0.3&0.4\\0.2&-0.2&-0.3\\-0.3&0.4&0.2\end{matrix}\right]\mathbf{y_{t-1}}+
\left[\begin{matrix}0.1&-0.2&0.5\\-0.4&0.1&-0.4\\0.1&0.3&-0.2\end{matrix}\right]\mathbf{y_{t-2}}+
\boldsymbol{\varepsilon_t}
}$$
$${\scriptsize
\boldsymbol{\varepsilon}_t\sim W.N.\left(\left[\begin{matrix}0.2^2&0&0\\0&0.2^2&0\\0&0&0.2^2\end{matrix}\right]\right)
}$$

とりあえずこのモデルに従うサンプルデータを作らないといけないですね。

まず、各係数行列や定数項などを変数として格納しておきます。

import numpy as np
import pandas as pd


# 各係数を変数に格納
A = np.array(
    [
        [1, 0, 0],
        [0.3, 1, 0],
        [-0.2, -0.3, 1]
    ]
)
A1 = np.array(
    [
        [-0.1, 0.3, 0.4],
        [0.2, -0.2, -0.3],
        [-0.3, 0.4, 0.2]
    ]
)
A2 = np.array(
    [
        [0.1, -0.2, 0.5],
        [-0.4, 0.1, -0.4],
        [0.1, 0.3, -0.2]
    ]
)
c = np.array([0.8, -0.5, 0.3])

これを使って、データを作ります。行列$A$の逆行列を両辺にかけて、モデルを誘導形に変形して順番に計算したらOKです。最初の2時点のデータは適当に設定して、初期のデータを捨ててます。

# 最初の2点のデータy_0, y_1を仮置き
df = pd.DataFrame(
    {
        "y0": [1, 1],
        "y1": [1, 1],
        "y2": [1, 1],
    }
)

# Aの逆行列
A_inv = np.linalg.inv(A)

# 乱数固定
np.random.seed(1)
for i in range(len(df), 550):
    df.loc[i] = A_inv@(c+A1@df.iloc[-1].T+A2@df.iloc[-2].T+np.random.normal(size=3)*0.2)

# 最初の方のデータを切り捨てる。
df = df.iloc[50:]
df.reset_index(inplace=True, drop=True)

これでデータができました。グラフは省略していますが、plotしてみると定常であることがわかります。

データができたので、statsmodelsで推定してみましょう。ドキュメントはこちらです。
参考: statsmodels.tsa.vector_ar.svar_model.SVAR — statsmodels
推定結果についてはこちら。
参考: statsmodels.tsa.vector_ar.svar_model.SVARResults — statsmodels

これの使い方は結構特殊です。左辺の行列のうち、推定したいパラメーターを文字列”E”とした行列を作って一緒に渡してあげる必要があります。対角成分は1です。前回の記事で書きましたが、各過程が外生性が高い順に並んでると仮定しているんので、下三角行列になるようにしています。具体的には次のようなコードになります。

import statsmodels.api as sm


# 上式のAの中で求めたい要素を"E"とした行列を生成しモデルに渡す。
A_param = np.array([[1, 0, 0], ["E", 1, 0], ["E", "E", 1]])

svar_model = sm.tsa.SVAR(df, svar_type="A", A=A_param)
svar_result = svar_model.fit(maxlags=2)

はい、これでできました。推定されたパラメーターが変数svar_resultに入っていますので、順番に見ていきましょう。

# 左辺の行列A
print(svar_result.A)
"""
[[ 1.          0.          0.        ]
 [ 0.34089688  1.          0.        ]
 [-0.20440697 -0.31127542  1.        ]]
"""

# 右辺の係数はcoefsにまとめて入っている
print(svar_result.coefs)
"""
[[[-0.12822445  0.30295997  0.43086499]
  [ 0.28516025 -0.30381663 -0.40460447]
  [-0.25660898  0.38056812  0.17305955]]

 [[ 0.05884724 -0.18589909  0.52278365]
  [-0.40304548  0.11406067 -0.60635189]
  [-0.07409054  0.28827388 -0.23204729]]]
"""

# 定数項
print(svar_result.intercept)
"""
[ 0.86592879 -0.82006429  0.28888304]
"""

そこそこの精度で推定できていますね。

VARモデルの時と同様に、summary()メソッドで結果を一覧表示することもできます。ただ、statsmodelsの現在のバージョン(0.13.5)のバグだと思うのですが、k_exog_userってプロパティを手動で設定しておかないと、AttributeError: ‘SVARResults’ object has no attribute ‘k_exog_user’ってエラーが出ます。とりあえず0か何か突っ込んで実行しましょう。

svar_result.k_exog_user=0
svar_result.summary()
"""
  Summary of Regression Results   
==================================
Model:                        SVAR
Method:                        OLS
Date:           Mon, 27, Feb, 2023
Time:                     00:29:40
--------------------------------------------------------------------
No. of Equations:         3.00000    BIC:                   -9.36310
Nobs:                     498.000    HQIC:                  -9.47097
Log likelihood:           276.728    FPE:                7.18704e-05
AIC:                     -9.54065    Det(Omega_mle):     6.89230e-05
--------------------------------------------------------------------
Results for equation y0
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.865929         0.038253           22.637           0.000
L1.y0        -0.128224         0.042999           -2.982           0.003
L1.y1         0.302960         0.037971            7.979           0.000
L1.y2         0.430865         0.041409           10.405           0.000
L2.y0         0.058847         0.043158            1.364           0.173
L2.y1        -0.185899         0.040139           -4.631           0.000
L2.y2         0.522784         0.035344           14.791           0.000
========================================================================

Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const        -0.820064         0.043394          -18.898           0.000
L1.y0         0.285160         0.048778            5.846           0.000
L1.y1        -0.303817         0.043074           -7.053           0.000
L1.y2        -0.404604         0.046974           -8.613           0.000
L2.y0        -0.403045         0.048958           -8.233           0.000
L2.y1         0.114061         0.045533            2.505           0.012
L2.y2        -0.606352         0.040094          -15.123           0.000
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.288883         0.041383            6.981           0.000
L1.y0        -0.256609         0.046517           -5.516           0.000
L1.y1         0.380568         0.041078            9.265           0.000
L1.y2         0.173060         0.044797            3.863           0.000
L2.y0        -0.074091         0.046688           -1.587           0.113
L2.y1         0.288274         0.043423            6.639           0.000
L2.y2        -0.232047         0.038236           -6.069           0.000
========================================================================

Correlation matrix of residuals
            y0        y1        y2
y0    1.000000 -0.300511  0.090861
y1   -0.300511  1.000000  0.269623
y2    0.090861  0.269623  1.000000
"""

正直、SVARモデルを使う時って、左辺の係数行列$A$が一番注目するところだと思うのですが、その情報が出てこないってところがイマイチですね。このsummary()はVARモデルほどは使わないと思いました。バグが放置されているのものそのせいかな?

最後の結果出力だけイマイチでしたが、必要な値は各属性を直接見れば取れますし、そこそ手軽に使えるモデルではあったのでVARモデルで力不足に感じることがあったらこれも思い出してみてください。

構造VARモデルの紹介

久しぶりに時系列の話です。以前、ベクトル自己回帰モデル(VAR)というのを紹介しました。
参考: ベクトル自己回帰モデル

これは要するに、時系列データのベクトルを、それより前の時点のベクトルで回帰する(線形和として表現する)ことによって説明しようっていうモデルでした。

これを実際に業務で使おうとすると、非常に厄介な問題が発生します。それは、同じ時点での値どうしの間にも関係があるということです。

例えば、何かのECサイトの分析をしてて、サイトの訪問者数、会員登録数、売上、の3つの時系列データがあったとした場合、VARの観点で言うと、過去の訪問者数、過去の会員登録数、過去の売上、からその日の各値を予測しようって言うのがVARモデルです。

そうなった時に、いやいや、会員登録数は「当日の訪問者数」の影響を受けるし、売上は「当日の会員登録数」の影響を受けるでしょうとなります。VARモデルではそう言う影響が考慮できません。ここを改善したのが構造VARモデル(Structural Vector Autoregressive Model)です。

時系列分析でいつも引き合いに出している、沖本先生の「経済・ファイナンスデータの計量時系列分析」では、4.6節(99〜101ページ)でさらっと紹介されています。2ページくらいです。

具体的に、3変数で最大ラグが2のSVARモデルを書き出すと下のような形になります。

$$\left\{\begin{align}y_{0,t}&=c_0&&& +\phi_{0,0,1}y_{0,t-1}+\phi_{0,1,1}y_{1,t-1}+\phi_{0,2,1}y_{2,t-1}+\phi_{0,0,2}y_{0,t-2}+\phi_{0,1,2}y_{1,t-2}+\phi_{0,2,2}y_{2,t-2}+\varepsilon_0\\
y_{1,t}&=c_1&-\phi_{1,0,0}y_{0,t}&&+\phi_{1,0,1}y_{0,t-1}+\phi_{1,1,1}y_{1,t-1}+\phi_{1,2,1}y_{2,t-1}+\phi_{1,0,2}y_{0,t-2}+\phi_{1,1,2}y_{1,t-2}+\phi_{1,2,2}y_{2,t-2}+\varepsilon_1\\
y_{2,t}&=c_2&-\phi_{2,0,0}y_{0,t}&-\phi_{2,1,0}y_{1,t}&+\phi_{2,0,1}y_{0,t-1}+\phi_{2,1,1}y_{1,t-1}+\phi_{2,2,1}y_{2,t-1}+\phi_{2,0,2}y_{0,t-2}+\phi_{2,1,2}y_{1,t-2}+\phi_{2,2,2}y_{2,t-2}+\varepsilon_2\\\end{align}\right.$$

後ろの方、ブログ幅からはみ出しましたね。見ての通りそこそこ巨大なモデルになります。
この、$y_{1,t}$の予測に$y_{0,t}$が使われていたり、$y_{2,t}$の予測に$y_{0,t},y_{1,t}$が使われているのが、VARモデルとの違いです。

逆に、$y_{0,t}$の説明には$y_{1,y}$や$y_{2,t}$は用いられてはいません。

これは変数が外生性が高い順に並んでいることを仮定しているためです。これを仮定せず、相互に影響し合うようなモデルも研究されてはいるようなのですが、実データから係数を推定するのが非常に難しくなるので、SVARモデルを使うときはこの仮定を置いておいた方が良い、と言うよりこれが仮定できる時に利用を検討した方が良いでしょう。

通常は、時刻$t$の時点の項を左辺に移行して行列表記するようです。(そのため、上の例でも移項を想定してマイナスつけときました。)

移行して、行列、ベクトルを記号に置き換えていくと、下のような式になります。$\mathbf{D}$は$n\times n$の対角行列とします。要するに撹乱項はそれぞれ相関を持たないとします。

$$
\mathbf{B}_0\mathbf{y}_t=\mathbf{c}+\mathbf{B}_1\mathbf{y}_{t-1}+\cdots+\mathbf{B}_p\mathbf{y}_{t-p}+\boldsymbol{\varepsilon}_t,\ \ \
\boldsymbol{\varepsilon}_t\sim W.N.(\mathbf{D})
$$

この形の式を構造形(Structual form)と呼びます。VARモデルとの違いは、左辺に$\mathbf{B}_0$が掛かってることですね。

この構造形ですが、実際にデータがあった時にこのまま係数を推定することが難しいと言う問題があります。そのため、両辺に$\mathbf{B}_0$をかけて次の形の式を考えます。

$$
\mathbf{y}_t=\mathbf{B}_0^{-1}\mathbf{c}+\mathbf{B}_0^{-1}\mathbf{B}_1\mathbf{y}_{t-1}+\cdots+\mathbf{B}_0^{-1}\mathbf{B}_p\mathbf{y}_{t-p}+\mathbf{B}_0^{-1}\boldsymbol{\varepsilon}_t
$$

この形を、誘導形(reduced form)と呼びます。撹乱項にも逆行列がかかってるので、誘導形の撹乱項の各成分には相関が生まれている点に気をつけてください。

この誘導形はVARモデルと全く同じなので、VARモデルを推定するのと同じ方法でパラメーターを推定することができます。

そして、ここからが問題なのですが、誘導形から構造形を求めることはそう簡単ではありません。(逆に構造形が事前にわかっていた場合に、誘導形を求めることは簡単です。上でやった通り、両辺に逆行列かけるだけだからです。)

その難しさの原因を説明をします。上の方で、変数が外生性が高い順に並んでると仮定して$\mathbf{B}$の一部の成分が0であることを仮定していましたが、もしその仮定がなかったとしましょう。すると、構造形の方が誘導形よりパラメーターが多くなってしまうのです。

構造形の方は、$\mathbf{y}_0$の係数の行列が、対角成分は1とわかってるので残り$n(n-1)$個、右辺は定数項が$n$個、右辺の各行列が全部で$pn^2$個、撹乱項が$n$個ので、合わせて、$n(n-1)+n+pn^2+n=n(n+1)+pn^2$個のパラメータを持ちます。

一方で誘導形の方は、左辺のパラメーターこそ消えますが、右辺は定数項が$n$個、右辺の各行列が全部で$pn^2$個、撹乱項が$n(n+1)/2$個で、合計$n+pn^2+n(n+1)/2$個しかパラメーターを持ちません。

その差、$n(n-1)/2$個だけ、構造形がパラメーターが多いので、誘導形が定まった時に構造形が決められなくなってしまいます。

そこで、誘導形から構造形を一意に定めるために、$n(n-1)/2$個の制約を課す必要が発生します。その制約として、$\mathbf{B}_0$の上三角成分(対角成分より上)を0と決めてしまうのが、変数が外生性が高い順に並んでいると言う仮定です。

個人的には、この仮定があったとしても誘導形から構造形を導くのってそこそこ難しいように感じますが、それでも理論上は算出が可能になりますね。

具体的にデータを用意してPythonを使ってSVARのパラメーター推定をやってと記事を続ける予定だったのですがここまでで結構長くなってしまったのと、数式のせいでそこそこ疲れてしまったので、それは次回の記事に回そうと思います。

VARの欠点をクリアしたモデルではありますが見ての通り巨大なので、かなりデータが多くないと推定がうまくいかなかったりして、なかなか期待ほど活躍しないのですが、こう言うモデルもあるってことを認識しておくとどこかで役に立つのかなと思います。

SciPyで数列の極大値や極小値を求める

時系列データを分析している中で、極大値や極小値を特定したいケースは稀にあります。

極大値/極小値というのは、要は局所的な最大値/最小値のことで、その値の周囲(前後)の値と比較して最大だったり最小だったりする要素のことです。(とても雑な説明。もっと正確な説明はWikipediaの極値のページを参照。)

例えば、[5, 4, 3, 2, 3, 4, 3, 2, 1] みたいな数列があった時、一番値が大きいのは先頭の5なので、これが最大値(極大値でもある)ですが、6番目の4もその近くだけ見ると、[2, 3, 4, 3, 2]となっていて前後の値より大きいので、この6番目の4が極大値です。

今の説明の通りのコードを書いて数列の前後の値と比較して判定したら極大値も極小値も見つかるのですが、それを一発でやってくれるメソッドがSciPyにあるよ、ってのが今回の記事です。

使うのはargrelmin/ argrelmax です。ドキュメントは以下。
scipy.signal.argrelmin — SciPy v1.10.0 Manual
scipy.signal.argrelmax — SciPy v1.10.0 Manual

minとmaxの違いは極小値か極大値かの違いなので、以下の説明はargrelmaxの方でやっていきますね。

ちょっと適当な数列を一個用意して実行してみます。(さっきの例に値を付け足して長くしたものです。)

import numpy as np
from scipy.signal import argrelmax

# サンプルデータを用意
data = np.array([5, 4, 3, 2, 3, 4, 3, 2, 1, 8, 1])

# 極大値のindex
print(argrelmax(data))
# (array([5, 9]),)

# 極大値の値
print(data[argrelmax(data)])
# [4 8]

indexが5(インデックスは0始まりなので6番目の要素)である4と、indexが9(同様に10番目の要素)である8が検出されました。

printした結果を見ていただくと分かる通り、argrelmaxはindexが入ったarrayを0番目の要素に持つタプル、という特殊な形で結果を返してくれます。慣れないとトリッキーに見えますが、それをそのまま使うと極値の値を取り出せるので便利です。

デフォルトでは、直前直後の値だけを見て極大値極小値が判定されますが、例えばノイズを含むデータなどでは実用上検出が多すぎることもあります。
その場合、order(デフォルト1、1以上の整数のみ指定可能)という引数を使うことで、前後何個の要素と比較するかを指定できます。

order=3 とすると、前の3個の値と、後ろの3個の値、合計6個より大きい場合に極大値として判定されます。

data = np.array([2, 1, 1, 4, 1, 1, 5, 1])

# index3の4 と index6の5が極大値として検出される。
print(argrelmax(data, order=1))
# (array([3, 6]),)

# order=3とすると、index6の5だけが極大値として検出される。
print(argrelmax(data, order=3))
# (array([6]),)

上記の例でもう一個着目して欲しい点があります。order=1の時、先頭の2は極大値として検出されませんでした。デフォルトでは、orderの値に関係なく前後に最低1個の要素がないと対象にならないようです。そして、order=3の場合も、後ろから2番目の5が検出されています。orderで指定した数に足りなくても前後に1個以上あれば良いようです。

この、端に関する挙動はmodeという引数で指定できます。デフォルトは”clip”で、これは両端の値は極値として扱われません。ここに”wrap”を指定すると、両端の値も対象になります。

# index=0の2も対象になった
print(argrelmax(data, order=1, mode="wrap"))
# (array([0, 3, 6]),)

もう一つ気をつけないといけないのは、ドキュメントに書かれている通り、前後の値より真に大きくないと極値として扱われません。以下のように前後の値と一致したらダメということです。

data = np.array([1, 1, 2, 3, 3, 2, 1])
print(argrelmax(data))
# (array([], dtype=int64),)

以上が1次元のデータに対する使い方になります。

さて、このargrelmin/ argrelmaxですが、2次元以上のデータに対しても使えます。ドキュメントには2次元の例が載っていますが、3次元でも4次元でもいけます。

2次元、要するに行列形式のデータに対して使ったら、上下左右、できれば斜めも考慮した8方向の隣接データと比較して極大値/極小値を出してくれるのかな?と期待したのですがそういう動きはしておらず、軸(axis)を1個固定してその軸に沿った1次元データとして取り出してそれぞれに対して極大値/極小値の検索をやるようです。方向はaxis引数(デフォルト0)で指定します。

ちょっとでたらめに作ったデータでやってみます。

data = np.array([
        [1, 4, 0, 9, 0, 3, 2],
        [4, 0, 0, 1, 2, 3, 7],
        [2, 9, 2, 0, 9, 0, 7],
        [0, 0, 7, 9, 6, 3, 1],
        [0, 4, 4, 7, 2, 8, 3]
    ])

# axis省略は0と同じ。
print(argrelmax(data))
# (array([1, 2, 2, 3, 3]), array([0, 1, 4, 2, 3]))

print(argrelmax(data, axis=1))
# (array([0, 0, 0, 2, 2, 3, 4, 4]), array([1, 3, 5, 1, 4, 3, 3, 5]))

結果の読み方が慣れないと分かりにくいですが、インデックスの1次元目の配列、2次元目の配列のタプルとして帰ってきてます。

要するに、 (array([1, 2, 2, 3, 3]), array([0, 1, 4, 2, 3])) というのは [1, 0], [2, 1], [2, 4], [3, 2], [3, 3]が極大値だった、という意味です。

そして、またこれも分かりにくいのですが、axis=0 の時、この行列の各列を取り出して極大値を探しています。[1, 0]の値は4 ですが、これはdata[:, 0] = [1, 4, 2, 0, 0] の極大値として検出されており、[3, 2]の7 は data[:, 2] = [0, 0, 2, 7, 4] の極大値として検出されています。

スライスした時に:(コロン)になる次元をaxis引数で指定していると考えたら良いでしょうか。

引数を省略してた時の挙動が想定と違うというか、axis=1を指定した時の方がデフォルトっぽい動きしていますね。こちらは、
(array([0, 0, 0, 2, 2, 3, 4, 4]), array([1, 3, 5, 1, 4, 3, 3, 5]))
が結果として帰りますが、こちらも同様に[0, 1], [0, 3], [0, 5], ….(略) が極大値として検出されています。そしてこれは data[0, :] = [1, 4, 0, 9, 0, 3, 2] の極大値です。

滅多に使わない関数ですしさらにこれを多次元データに使うというのも稀だと思うので、完璧に理解し記憶して使いこなすというよりも、必要になった時に挙動をテストしながら使うのが現実的ではないでしょうか。

M2搭載のMacBookにPython環境構築 (2023年02月時点)

最近、私物のMacBookを買い替えました。買ったのは 「MacBook Air M2 2022」です。(以前はPro使ってましたが、開発環境はAWS上にもあるのでローカルはAirで十分かなと思って変えました。)

さて最近のMacBookは、CPUがIntel製ではなく、Apple製のものになっています。この影響で特にM1登場当初は環境構築で苦労された人もたくさんいたようですし、すでに多くの記事が書かれてナレッジがシェアされています。僕自身、それが原因で調子が悪かった先代のMBPをなかなか買い替えなかったというのもあります。

しかし現在ではOSSコミュニティーの皆さんの尽力のおかげで環境はどんどん改善しており、僕自身はそこまで大きな苦労なくPython環境を構築できました。とはいえ、あくまでも私用端末なので、業務で使ってる端末に比べると入れたライブラリ等がかなり少ないのですが。それでも、初期の頃つまづいた報告が多かったnumpyやmatplotlibなどは入れれることを確認しています。

現時点ではこんな感じですよ、ってことで誰かの参考になると思うので手順を書いていきます。

前提:
– シェルはzsh
– pyenv利用
– Anaconda/ Minicondaは使わない
– Python 3.11.1 をインストール

では、やっていきます。

1. Homebrewインストール

pyenvを導入するために、まずHomebrewを入れます。
公式サイトに飛んでインストールコマンドを実行します。一応コマンドはこのブログにも載せますが、昔とコマンドが変わっていますし、今後変わる可能性もありますし実行する時にちゃんと公式サイトを確認したほうがいいですね。

公式サイト: macOS(またはLinux)用パッケージマネージャー — Homebrew

# サイト掲載のインストールコマンド
$ /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

以前だったらこれを実行して、メッセージに従って何度かENTERキー押して進めていけば終わりだったのですが、今は最後に以下のようなメッセージが表示されます。(正直、見落としがちなので公式サイトにも書いておいて欲しかった。)

Warning: /opt/homebrew/bin is not in your PATH.
  Instructions on how to configure your shell for Homebrew
  can be found in the 'Next steps' section below.
==> Installation successful!
# -- 中略 --
==> Next steps:
- Run these three commands in your terminal to add Homebrew to your PATH:
    echo '# Set PATH, MANPATH, etc., for Homebrew.' >> /Users/{ユーザー名}/.zprofile
    echo 'eval "$(/opt/homebrew/bin/brew shellenv)"' >> /Users/{ユーザー名}/.zprofile
    eval "$(/opt/homebrew/bin/brew shellenv)"
- Run brew help to get started
- Further documentation:
    https://docs.brew.sh

要するにPathが通ってないから通して、とのことです。

Next stepsにある3行を実行します。ただ、1行目はただのコメント文を残すだけですし、3行目のeval はターミナル/シェルを再実行すれば良いだけなので、必須なのは2行目だけです。
この記事ではマスクしてますが、{ユーザー名}部分、ちゃんと端末のユーザー名がメッセージに表示されてます。気が利きますね。

$ echo '# Set PATH, MANPATH, etc., for Homebrew.' >> /Users/{ユーザー名}/.zprofile
$ echo 'eval "$(/opt/homebrew/bin/brew shellenv)"' >> /Users/{ユーザー名}/.zprofile
$ eval "$(/opt/homebrew/bin/brew shellenv)"

ちなみに、これをやると環境変数PATHにhobebrew関係のパスが追加されます。気になる人は実行前後で見比べてみましょう。

2. pyenvインストール

Homebrewが入ったら次はpyenvです。ドキュメントにHomebrewを使って入れる専用のセクションがあるのでそれに従います。

ドキュメント: https://github.com/pyenv/pyenv/blob/master/README.md#homebrew-in-macos

$ brew update
$ brew install pyenv

続いて、シェルの設定です。こっちにあります。zshの方をやります。
参考: https://github.com/pyenv/pyenv/blob/master/README.md#set-up-your-shell-environment-for-pyenv

$ echo 'export PYENV_ROOT="$HOME/.pyenv"' >> ~/.zshrc
$ echo 'command -v pyenv >/dev/null || export PATH="$PYENV_ROOT/bin:$PATH"' >> ~/.zshrc
$ echo 'eval "$(pyenv init -)"' >> ~/.zshrc

3. Pythonインストール

pyenvが入ったらPythonを入れます。ちゃちゃっと入れたいところですが、そのまいれたら、
WARNING: The Python lzma extension was not compiled. Missing the lzma lib?
という警告が出ました。xzってのを先に入れておくとこれが出ないのでやっておきましょう。homebrewで入ります。

$ brew install xz

ここまでできたらPythonインストールとバージョン切り替えです。

入れたいバージョンがはっきり決まっているならそれをそのまま入れたら良いのですが、僕は毎回その時点でインストール可能なバージョンの一覧を眺めてから決めています。今回は3.11.1を入れました。

# インストール可能なバージョン一覧表示 
$ pyenv install -l 
# インストール 
$ pyenv install 3.11.1
# バージョン切り替え
$ pyenv global 3.11.1

そして、ライブラリを入れていきます。requirements.txt作って一括で入れてもいいのですが、ちょっと怖かったので1個ずつ恐る恐る入れていきましたが、概ねすんなり入っていくようです。正直、Pythonだけで書かれてるライブラリたちは比較的安心なのですが、numpy等のC言語が使われているライブラリは不安でしたね。

$ pip install jupyterlab
$ pip install numpy
$ pip install pandas
$ pip install scipy
$ pip install matplotlib
$ pip install scikit-learn
$ pip install lxml
$ pip install requests
$ pip install beautifulsoup4
# 以下略

以上のようにして、M2 MacでもPythonが使えるようになりました。

対応を進めていただいた開発者の皆様、ありがとうございました。

LINE NotifyでPythonから自分にメッセージを送る

前回に続いて通知を作る話です。
今回はLINE通知を作ります。

この記事は、企業が運用している本格的な公式アカウントやBotサービスのようなものではなく、個人的に運用しているサーバーのバッチでエラーが起きたときなどに自分宛に通知するという小規模な利用を想定して書きます。

利用するサービスはこちらの LINE Notify です。
参考: https://notify-bot.line.me/ja/

この種のSNSに付随したサービスに対しては認証周りで面倒なコードを書かないといけないイメージがあったのですが、LINE Notifyは非常にコンパクトな実装で手軽に使えました。

事前準備として、こちらのサービスからトークンを入手します。

これ専用のアカウントは必要なく、普段使っているLINEのアカウントでログインできます。右上のログインリンクからログインしましょう。(LINE認証用のコードが表示され、LINEアプリから入力が必要なので、スマホも用意しておきましょう。)

ログインしたら、右上のログインリンクが自分のLINEアカウント名になっているので、それを押してマイページへ遷移します。
そして、「アクセストークンの発行(開発者向け) 」のところから「トークンを発行」ボタンをクリックします。

トークン名を入力し、通知を送るトークルームを選択して、「発酵する」ボタンをクリックするとトークンが発行されて1度だけ表示されます。もう二度と表示されないのでこの時点で確実に記録しておきましょう。

先に書いておきますが、通知を実装した後実際にLINEに届くメッセージは、
[トークン名] メッセージ本文
というフォーマットになります。トークン名が長いと毎回邪魔なのでコンパクトでわかりやすい名前にしておきましょう。

さて、トークンが発行できたらこれを使ってみます。
ドキュメントはこちらです。
参考: LINE Notify API Document
このドキュメントの「通知系」のところにある、https://notify-api.line.me/api/notify が通知を送るAPIです。

リクエストパラメーターがいろいろ書かれていますが、「必須」と指定されているのはmessageだけなので非常に簡単に使えます。

CURLで動かすサンプルコードもあるのでちょっとやってみましょう。{自分のトークン}の部分は先ほど発行したトークンを入れてください。

# コマンド
$ curl -X POST -H 'Authorization: Bearer {自分のトークン}' -F 'message=CURLで通知' https://notify-api.line.me/api/notify

# 結果
{"status":200,"message":"ok"}

これで、「[トークン名] CURLで通知」というメッセージが、 LINE Notify のアカウントから届きます。

あとは、このcurlコマンドをPythonに書き直していきましょう。使うライブラリはrequestsあたりで良いと思います。

import requests


line_notify_token = "{自分のトークン}"
api_url = "https://notify-api.line.me/api/notify"
message = "メッセージ本文"
headers = {"Authorization": f"Bearer {line_notify_token}"}
data = {"message": message}
requests.post(
    api_url,
    headers = headers,
    data = data,
)

たったこれだけで、LINEにメッセージが届きます。

もう少し丁寧に実装するなら、postの戻り値のstatusコードを確認してエラー処理を入れたりするとよさそうですね。

LINEのインターフェース的に、あまりにも長文を送ったりするのには適さず、長文になるなら先日のメール通知の方が良いかなと思うのですが、
メールよりLINEの方が通知に気付きやすいので、速報性が必要な場面で重宝しそうです。

PythonでYahooメールのアカウントからメール送信

個人的に運用しているソフトウェアに通知機能を作りたかったので、メールの送信方法を調べました。本当はSlack通知とかの方が使いやすいのですが、最近のSlackの無料プランは一定期間でメッセージが消えるなどイケてないですからね。

メールの利用を検討し出した当初はAmazon SESを使おうかとも思っていたのですが、標準ライブラリだけで実装できることと、料金もかからないのでPythonでやることにしました。

使用するライブラリは以下の二つです。
– SMTPプロトコルクライアント – smtplib
– メールメッセージの管理 – email
email の方は 使用例のページを見た方がいいです。

使い方はYahooメールやGmail, Outlookなどアカウントによって微妙に違うので今回はYahooメールを例に取り上げて説明します。

メールサーバーやポート番号の情報はこちらにあります。
参考: メールソフトで送受信するには(Yahoo!メールアドレス、@ymail.ne.jpアドレスの場合)

必要な情報は以下の2つです。
– 送信メール(SMTP)サーバー : smtp.mail.yahoo.co.jp
– 送信メール(SMTP)ポート番号 : 465

実はこの情報は、2020年8月に変わっていて、他所の古い技術記事等ではポート番号が違っていたりします。他サイトのコードをコピペしたが動かなかったという人は以下のアナウンスを読んでください。SMTPの通信方法がSSLになっているというのもライブラリで呼び出す関数が変わるので重要な点です。
参考: Yahoo!メールをより安全にご利用いただくためのメールソフト設定(送受信認証方式)変更のお願い

メールソフトで送受信するにはのページに記載がありますが、「Yahoo! JAPAN公式サービス以外からのアクセスも有効にする」の設定をやっておかないと動かないので気をつけてください。(とはいえ、普段スマホでメールを見れるようにしているのであれば設定済みだと思います。)

それではやっていきましょう。自分のスクリプトから自分のメールアドレス宛の通知としての利用を想定しているので飾りも何もないテキストメールを送ります。

まず、各変数に必要な値を格納しておきます。悪用は厳禁ですが、toだけでなくfromのメールアドレスも自由に設定できます。まともに使うならfromは自分のアドレスでしょう。

# Yahooのログイン情報
username = "{YahooのユーザーID}"
password = "{Yahooのパスワード}"

# メールアドレス情報
from_address = "{差出人のアドレス}@yahoo.co.jp"
to_address = "{宛先のアドレス}@yahoo.co.jp"

# SMTPサーバーの情報。値はYahooメールのヘルプページから取得。
smtp_host = "smtp.mail.yahoo.co.jp"
smtp_port = 465

続いて、メールの情報を作ります。smtpのドキュメント末尾では文字列で直接データを作ってますが、その直下の注釈でemailパッケージを推奨されているのでそちらに従います。

from email.message import EmailMessage


msg = EmailMessage()
msg.set_content("メール本文")
msg["Subject"] = "メールタイトル"
msg["From"] = from_address
msg["To"] = to_address

メールデータできたので、これを送信します。ドキュメントの一番下のサンプルコードではローカルのSMTPサーバーでメール送信しているのでログインも何もしていませんが、Yahooメールを使うなら最初にログインが必要です。smtplib.SMTP ではなく、smtplib.SMTP_SSLを使うのもポイントですね。コードは以下のようになります。

import smtplib


server = smtplib.SMTP_SSL(smtp_host, smtp_port)
server.login(username, password)
server.send_message(msg)
server.quit()

たったこれだけでメール送信が実装できました。

WordPressのログインページのURLを変更する

前回の記事で書いてる通り、機械的にログインを試みる攻撃を受けていたので対策を施しました。特定IPアドレスからの攻撃だったのでそのIPをブロックしようかと思ったのですが、他のIPに変えられるたびにやるのも面倒なので、Wordpressのセキュリティ策としてよく挙げられているログインURLの変更を実施しました。

初期設定のログインURLはWordpressのドキュメントを見れば分かっちゃいますからね。

方法はいろいろありますが、手軽な方法としてプラグインを使うことにしました。

選んだのは、 All In One WP Security です。もっとシンプルな、URLの変更に特化したやつとかもあるのですが今後別の対策を考える時があったら使いまわせるのがいいと思ったのでこれを選びました。(結果的にURL変更だけで攻撃が収まったのですが、対策前の時点ではそれだけで完了するかどうわかりませんでしたし。)

WordPressの管理画面からプラグインの画面を開き、新規追加から検索してインストールします。そして有効化します。

有効化したら左ペインのメニューに「WPセキュリテイ」というのができるのでここから設定します。

「総当たり攻撃」というカテゴリ(英語だとBrute force)の中の、Rename login page タブがログインURLの変更です。それを開きます。

Enable rename login page feature: のチェックボックスにチェックを入れ、Login page URL: のテキストボックスの中にこれから使うURLを設定しSaveするとそちらが新しいログインページになります。

これで、デフォルトのログインページにアクセスしてみてフォームが出てこないことを確認したら完成です。

WordPressは未ログイン状態でログイン後の管理画面のURLにアクセスするとログインURLにリダイレクトされたりするのですが、このツールでURLを変更するとその動作もなくなり、リダイレクトによってURLがバレるということも防いでくれています。気がききますね。

そして、肝心の攻撃に対する効果ですが、apacheのログを見ると該当の攻撃者に404エラーが出た段階でピタリとアクセスが止まってるのを確認できました。Lightsailのリソースも回復しており良い感じです。

今回の事象への対応はこれで完了ですが、時々不自然にアクセスが集中している時間があったりなどの不穏な動きはまぁまぁあるので必要に応じて今回導入したAll-In-One Securityの各機能を活用して対策して行こうと思います。