PythonでSQLite3を使う

前回の記事に引き続き、SQLite3の話です。

今回はPythonのコードの中でSQLite3を使っていきます。Pythonでは標準ライブラリにSQLite3用のインターフェースがあり、特に何か追加でインストールしなくてもSQLite3を使えます。MySQL等と比較したときの大きなアドバンテージですね。
参考: sqlite3 — DB-API 2.0 interface for SQLite databases — Python 3.11.3 documentation

使い方はドキュメントにあるので早速使っていきましょう。

DB(SQLiter3のDBはただのファイル)は前回の記事で作ったのがあるのでそれに接続します。存在しなかったら勝手に生成されるのでここはあまり難しいところではありません。

前回の記事を少しおさらいしておくと、DBファイル名は、mydatabase.dbで、userってテーブルができています。これの全レコードをSELECTするコードは次の様になります。

import sqlite3


# データベースに接続
con = sqlite3.connect('mydatabase.db')
cursor = con.cursor()
# SQLの実行と結果の表示
cursor.execute("SELECT * FROM user")
row = cursor.fetchall()
print(row)
# 以下結果
# [(1, 'Alice', 20), (2, 'Bob', 30), (3, 'Charlie', 40)]

簡単ですね。

結果はタプルのリストで得られました。列名(SELECT文ではアスタリスク指定にしている)の一覧が欲しいと思われますが、それは、coursor.descriptionという属性に入ってます。ただ、入り方がちょっとクセがあって、DB API の仕様との整合性のために、不要なNone が入ってます。

# 列名 + 6個の None という値が入ってる
print(cursor.description)
"""
(
    ('id', None, None, None, None, None, None),
    ('name', None, None, None, None, None, None),
    ('age', None, None, None, None, None, None)
)
"""

# 列名の一覧を得るコードの例
print([row[0] for row in cursor.description])
# ['id', 'name', 'age']

これだとちょっと使いにくいということもあるともうのですが、row_factory というのを使うと、列名: 値 の辞書で結果を得ることができます。ドキュメントに「How to create and use row factories」ってしょうがあるのでそれを参考にやってみます。

con = sqlite3.connect('mydatabase.db')
# row_factoryを指定
con.row_factory = sqlite3.Row
cursor = con.cursor()
cursor.execute("SELECT * FROM user")
row = cursor.fetchall()
print(row)
# [<sqlite3.Row object at 0x109b4de70>, <sqlite3.Row object at 0x109b4dae0>, <sqlite3.Row object at 0x109b4fc70>]

# これらはそれぞれ辞書形でデータが入っているので次の様にして値を取り出せる。
print(row[0]["name"])
# Alice
print(row[0]["age"])
# 20


# 上記の状態だと、DataFrameへの変換が簡単
import pandas as pd


print(pd.DataFrame(row))
"""
   0        1   2
0  1    Alice  20
1  2      Bob  30
2  3  Charlie  40
"""

これでSELECTはできましたね。

次はINSERT等のデータの編集の話です。INSERTだけ取り上げますが、UPDATEもDELETEも同様です。

基本的にexecuteでクエリを実行できるのですが、端末でSQLite3を使う時と違って、明示的にコミットしないと、セッションを切って繋ぎ直したらデータが消えてしまいます。commit()を忘れない様にしましょう。

con = sqlite3.connect('mydatabase.db')
cursor = con.cursor()
# 1行挿入
cursor.execute("INSERT INTO user (name, age) VALUES (?, ?)", ["Daniel", 35])
# コミット
con.commit()
# レコードが1レコード増えてる
cursor.execute("SELECT * FROM user")
row = cursor.fetchall()
print(row)
# [(1, 'Alice', 20), (2, 'Bob', 30), (3, 'Charlie', 40), (4, 'Daniel', 35)]

この他 CREATE TABLE等々も実行できますので以上の知識で通常の使用はできると思います。

ただ、Pythonで使う場合と端末で直接使う場合で異なる点もあります。それが、.(ドット)で始まるSQLite3の固有コマンドたちが使えないってことです。

.table が使えないのでテーブルの一覧を得るには代わりに次の様にします。

# テーブル名の一覧を取得するSQL文を実行
cursor.execute("SELECT name FROM sqlite_master")
print(cursor.fetchall())
# [('user',)]

そして、特定のテーブルのスキーマを確認するクエリはこちらです。userテーブルしか作ってないので、userテーブルで実験します。

# スキーマを取得するSQL
cursor.execute("SELECT sql FROM sqlite_master WHERE name = 'user'")
print(cursor.fetchone()[0])
"""
CREATE TABLE user (
    id INTEGER PRIMARY KEY,
    name TEXT,
    age INTEGER
)
"""

どちらも sqlite_master ってテーブルからデータを引っ張ってきてることからわかる通り、メタデータ的なテーブルが存在してその中に情報があるということです。

ちなみに、この sqlite_master に含まれる列名の一覧は次の様になってます。 sqlite_master の中に sqlite_master の情報はなさそうなので、*で全列SELECTして列名の一覧を取得しています。

cursor.execute("SELECT * FROM sqlite_master WHERE name = 'user'")
print([row[0] for row in cursor.description])
# ['type', 'name', 'tbl_name', 'rootpage', 'sql']

以上で、PythonでSQLite3を使える様になると思います。

SQLite3入門

SQLite3の概要

SQLite3というのは軽量なRDBMSです。C言語で書かれてます。特徴としては、データベースをファイルとして扱い、サーバーが不要って点が挙げられます。

リレーショナルデータベースをちょっと使いたいけど、わざわざMySQLを導入するほどではないとか既存のDB環境は汚したくないってときに重宝します。

ただ、当然ですがMySQLではないので、SQL以外のコマンド群、MySQLで言えばSHOWで始まる各種コマンドなどが無く、その代わりにSQLite3専用のコマンドが用意されていたりして、最初は少し戸惑ったので記事としてまとめておきます。今回はとりあえず端末で直接使う場面を想定して書きます。次回の記事で、Pythonライブラリから使う方法を書く予定です。

インストール方法

macの場合。最初からインストールされています。
バージョンを上げたい場合は、Homebrewでもインストールできるようです。

$ brew install sqlite3

Linux (EC2のAmazon Linux2 を想定) の場合、こちらも最初から入っていますが、yumで入れることもできます。

$ sudo yum update
$ sudo yum install sqlite

ディストリビューションによってはyumでは無く、apt-get等別のコマンドのことがあるので確認して入れたほうが良いでしょう。

Windowsの場合、公式サイトのダウンロードページでバイナリファイルが配布されているのでダンロードして使います。

SQLite3の起動と終了

MySQLとかだと接続情報(ユーザー名やパスワード、外部であればホスト名など)を含む長いコマンドを打ち込んで接続しますが、SQLite3はDBファイル名を指定して実行するだけでそのファイルに接続してDBとして使えるようになります。1ファイルが1データベースで、その中にテーブルを作っていくイメージです。MySQLは一つのRDBMSに複数のデータベース(スキーマ)を作れるのでこの点が違いますね。

mydatabase.db というファイルで起動する場合のコマンドは下記です。既存ファイルを指定しても良いし、無ければ新規に作成されます。
sqlite3 を起動するとプロンプトが sqlite> に変わります。

$ sqlite3 mydatabase.db
SQLite version 3.39.5 2022-10-14 20:58:05
Enter ".help" for usage hints.
sqlite>

終了するには、 .exit か .quit を入力します。 .exitの方はリターンコードを一緒に指定することもできます。先頭に.(ドット)が必要なので注意してください。

sqlite> -- 終了コマンド2種類。どちらかで閉じる。
sqlite> .quit
sqlite> .exit

個人的なぼやきですが、sqlite3の入門記事を書く人はSELECT文とかより先に終了方法を書くべきだと思っています。MySQLのノリでquitとかexitとか入れると syntax error が起きます。vimの終了方法並みに重要な知識です。

SQLite3におけるSQLの実行

SQLの実行は特に変わったことがありません。MySQLとほぼ同じように実行できます。

テーブルを作って、データを入れて、それを選択してみましょう。

-- テーブルの作成
sqlite> CREATE TABLE user (
    id INTEGER PRIMARY KEY,
    name TEXT,
    age INTEGER
);

-- データの挿入
sqlite> INSERT INTO user (name, age) VALUES ('Alice', 20);
sqlite> INSERT INTO user (name, age) VALUES ('Bob', 30);
sqlite> INSERT INTO user (name, age) VALUES ('Charlie', 40);

-- データの検索
sqlite> SELECT * FROM user;
-- 以下結果
1|Alice|20
2|Bob|30
3|Charlie|40

データ挿入時に id は指定しませんでしたが勝手に連番入れてくれてます。これはsqlite3の特徴の一つで、すべてのテーブルにはrowidっていう一意の行idが振られていて、 INTEGER PRIMARY KEY の列を作るとその列はrowid列の別名となり、自動的に振られたidが格納されるようになるのです。

DELETE や UPDATEなども他のRDBMSと同様の構文で行えるのですがこの記事では省略します。

SQLite3のコマンド

冒頭に少し書いた通り、SQLite3では、SQL以外にも独自のコマンドがあります。その代わり他のRDBMS固有のコマンドは使えません。SHOWなんとか、と打ちたくなって動かないことに気づいたら思い出しましょう。

そのコマンド群は . (ドット) 始まりになっています。代表的なものをいつくか挙げておきます。

  • .help: コマンドの一覧と使用方法を表示
  • .tables: データベース内のテーブル一覧を表示
  • .schema: テーブルのスキーマを表示 (指定したテーブルのCREATE TABLE文が見れる)
  • .backup: データベースをバックアップする
  • .import: CSVファイルなどの外部ファイルからデータを読み込み、テーブルに挿入する
  • .quitまたは.exit: SQLite3コマンドラインツールを終了する
  • .headers on|off: Selectした結果のヘッダーとして列名を表示するかどうか切り替える。初期設定はoff

一番最初に書いた .help で コマンドを一覧表示できるので一度一通りみておくと良いでしょう。(helpの起動方法が独自コマンドってこれもなかなか初心者に不親切な設計ですよ。)

まとめ

以上の内容を把握しておけば、最低限の用途でSQLite3を利用できるではと思います。

そして、最低限の用途で済まなくなったらSQLite3について詳しくなるよりも他のRDBMSへの移行を検討していくのが良いかなと考えています。簡易版RDBMSですからね。

matplotlibでオリジナルのカラーマップを作る

matplotlibでオリジナルの配色を使いたかったのでその方法を調べました。基本的にグラフを書くときにそれぞれのメソッドの引数で逐一RGBの値を指定していっても独自カラーを使うことはできるのですが、カラーマップとして作成すると、既存の配色と同じような雰囲気で使うことができるので便利です。

既存のカラーマップでも利用されている、以下の二つのクラスを使って作成できます。
matplotlib.colors.ListedColormap (N色の色を用意して選択して使う)
matplotlib.colors.LinearSegmentedColormap (線形に連続的に変化する色を使う)

上記のリンクがそれぞれのクラスのドキュメントですが、別途カラーマップ作成のページもあります。
Creating Colormaps in Matplotlib — Matplotlib 3.7.1 documentation

まず、ListedColormapのほうから試してみましょう。これはもう非常に単純で、色のリストを渡してインスタンスを作れば完成です。また、name引数で名前をつけておけます。省略するとname=’from_list’となり、わかりにくいので何かつけましょう。名前をつけておくとregisterってメソッドでカラーマップの一覧に登録することができ、ライブラリの既存のカラーマップと同じように名前で呼び出せるようになります。

それでは適当に、赤、緑、青、黄色、金、銀の6色で作ってみます。色の指定はRGBコードの文字列、0~1の実数値3個のタプル、matplotlibで定義されている色であれば色の名前で定義できます。例として色々試すためにそれぞれやってみます。

from matplotlib.colors import ListedColormap
import matplotlib as mpl
import numpy as np

colors = [
    "#FF0000",  # 赤
    (0, 1, 0),  # 緑
    (0, 0, 1, 1),  # 青(不透明度付き)
    "yellow",
    "gold",
    "silver",
]
pokemon_cmap = ListedColormap(colors, name="pokemon")

# 配色リストに登録するとnameで使えるようになる。(任意)
mpl.colormaps.register(pokemon_cmap)

これで新しいカラーマップができました。試しに使ってみましょう。文字列でつけた名前を渡してもいいし、先ほど作ったカラーマップのインスタンスを渡しても良いです。

# データを作成
np.random.seed(0)
data = np.random.randn(10, 10)

# 作図
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1)
psm = ax.pcolormesh(data, cmap="pokemon")
# psm = ax.pcolormesh(data, cmap=pokemon_cmap)  # 名前でななくカラーマップインスタンスを渡す場合
fig.colorbar(psm, ax=ax)
plt.show()

このような図ができます。色が適当なので見やすくはありませんが、確かに先ほど指定した色たちが使われていますね。

次は連続的に値が変化するLinearSegmentedColormapです。これは色の指定がものすごくややこしいですね。segmentdataっていう、”red”, “green”, “blue”というRGBの3色のキーを持つ辞書データで、それぞのチャネルの値をどのように変化させていくのかを指定することになります。サンプルで作られているデータがこれ。

cdict = {'red':   [(0.0,  0.0, 0.0),
                   (0.5,  1.0, 1.0),
                   (1.0,  1.0, 1.0)],

         'green': [(0.0,  0.0, 0.0),
                   (0.25, 0.0, 0.0),
                   (0.75, 1.0, 1.0),
                   (1.0,  1.0, 1.0)],

         'blue':  [(0.0,  0.0, 0.0),
                   (0.5,  0.0, 0.0),
                   (1.0,  1.0, 1.0)]}

各キーの値ですが3値のタプルのリストになっています。

リスト内のi番目のタプルの値を(x[i], y0[i], y1[i]) とすると、x[i]は0から1までの単調に増加する数列になっている必要があります。上の例で言えば、x[0]=0でx[1]=0.5でx[2]=1なので、0~0.5と0.5~1の範囲のredチャンネルの値の変化がy0,y1で定義されています。

値がx[i]からx[i+1]に変化するに当たって、そのチャンネルの値がy1[i]からy0[i+1]へと変化します。上記の”red”の例で言えば値が0〜0.5に推移するに合わせて0から1(255段階で言えば
255)へと変化し、値が0.5~1と推移する間はずっと1です。

y0[i]とy1[i]に異なる値をセットするとx[i]のところで不連続な変化も起こせるみたいですね。

あと、連続値とは言え実装上はN段階の色になる(そのNが大きく、デフォルトで256段階なので滑らかに変化してるように見える)ので、そのNも指定できます。

とりあえず上のcdictデータでやってみましょう。

newcmp = LinearSegmentedColormap('testCmap', segmentdata=cdict, N=256)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

pms = ax.pcolormesh(data, cmap=newcmp)
fig.colorbar(pms, ax=ax)
plt.show()

結果がこちら。

値が小さい時はRGBが全部0なので黒くて、徐々に赤みを増し、途中から緑の値も増え始めるので黄色に向かって、さらに途中から青も増え始めるので最大値付近では白になっていってます。わかりやすいとは言い難いですがなんとなく挙動が理解できてきました。

流石にこのやり方で細かい挙動を作っていくのは大変なので、LinearSegmentedColormap.from_listっていうstaticメソッドがあり、こちらを使うとただ色のリスト渡すだけで等間隔にカラーマップ作ってくれるようです。

また、その時zipを使って区間の調整もできます。使用例の図は省略しますが次のように書きます。(ドキュメントからそのまま引用。)

colors = ["darkorange", "gold", "lawngreen", "lightseagreen"]
cmap1 = LinearSegmentedColormap.from_list("mycmap", colors)

nodes = [0.0, 0.4, 0.8, 1.0]
cmap2 = LinearSegmentedColormap.from_list("mycmap", list(zip(nodes, colors)))

これのほうがずっと楽ですね。

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が使えるようになりました。

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