端末内の各フォルダに散らばっているファイルやドットファイルをGitで管理する

今回の記事は技術的には特に難しい話はなく、Gitのちょっとした応用なもので本職のエンジニアの人たちの間ではもしかしたら常識なのかもしれませんが、自分にとっては革新的だったので紹介しておきます。

プロジェクトのコードやファイルを管理するためにGitは普通に使っていると思いますが、何か大きなプロジェクトではなくて、端末(Macを想定)内の各ディレクトリに散っているファイルを個別にバージョン管理したくなることってないですか。

例えば、ホームディレクトリにある.vimrcなど.(ドット)始まりの隠しファイルとか、自分がよく使うコードをまとめた自作モジュールとか、頻繁に書くSQLの部品をまとめたメモとか、Pythonの環境構築に使う、requirements.txtなどもそうですね。
このファイルはGitで管理したいけど、それぞれの配置場所を個別にリポジトリにするのは面倒だし、同じディレクトリ内にGit管理が適さない属性のファイルもたくさんあるなぁ、っていう状況です。

このように、配置場所が散らばったファイルを1個のリポジトリで管理する方法があることを最近知りました。

やり方は簡単で、どこかに一つだけリポジトリを作り、その配下に各所に散らばっているファイルを集めてそれをgit管理し、元のフォルダにはシンボリックリンクを貼ると良いです。

まぁ、普通にディレクトリをどこかに掘って、リポジトリを作ります。
最初のコミットは空コミットにしておきましょう。

% mkdir my_files
% cd my_files
% git init
% git commit -m "first commit" --allow-empty

そして、このディレクトリへgit管理したいファイルたちを集めて、元ディレクトリへシンボリックリンク貼ってきます。
参考: ハードリンクとシンボリックリンク

% mv {元のファイルパス} {リポジトリ内のファイルパス}
% ln -s {リポジトリ内のファイルパス} {元のファイルパス}

どのファイルをどこにリンクしているかは、それはそれでREADME.md ファイルかどこかに一緒に記録しておくと良いでしょう。

ちなみに、リポジトリのルート直下に全ファイルまとめておくのはお勧めしません。dotファイルのディレクトリとかPythonモジュールのディレクトリとか切ってリポジトリ内を適切に整理しましょう。

こうすると、一つのリポジトリに各所に置いてあったファイルの実体が集まるのでgitでまとめて管理できるようになります。

そして本来の配置場所にはシンボリックリンクが貼られているので今までと全く同じように使用することができます。

注意点としては、ドットファイルの中でも特にセキュアな認証情報などを環境変数に設定するファイルを管理する場合の流出リスクですね。セキュアな情報はこの管理の対象外にするか、対象にするのであればgithub等外部のリポジトリへは上げない方が良いでしょう。(誤って公開リポジトリに上げてしまうと大事故に繋がり得ます)

J-Quants APIのページング処理に対応する

久々にJ-Quants API の記事です。もう結構前の話(2023/06/16)の話ですが、J-Quants APIはデータ量の増加位に対応するためにページング処理というものが導入されました。
参考: お知らせ – J-Quants API の 過去のお知らせ部分見てください。

要するにAPIから取得できるデータの量が多い時に、全部のデータを一度では取得できず、一部分だけ取得できるって話ですね。

こちらについて利用方法を記事にしておきます。

ページング処理対応方法

詳しくはこちらをご参照ください。
参考: API共通の留意事項 – J-Quants API

レスポンスが帰ってきた時、結果にpagination_key が含まれていたらページング(ページネーション)が発生しており、そこで得られた結果は取得したかったデータの全量ではありません。
得られたpagination_keyの値を付与して再度リクエストすることで以降のデータを得ることができます。

サンプルコード参照してやってみましょう。
ちなみに、認証にidトークンが必要ですがその取得方法は僕の過去記事参照してください。
参考: J-Quants API の基本的な使い方
以下の記事では、 id_token って変数にすでにトークンが取得できているものとします。

import json
import requests
import pandas as pd


print(len(id_token))  # id_tokenは過去記事の方法ですでに取得してるとします。(文字数確認)
# 1107

# 特定の日付の4本値を取得する
date = "2024-03-15"
daily_quotes_url = f"https://api.jquants.com/v1/prices/daily_quotes?date={date}"
headers = {"Authorization": f"Bearer {id_token}"}
daily_quotes_result = requests.get(daily_quotes_url, headers=headers)

# レスポンスに、pagination_key が含まれていることが確認できる。
print(daily_quotes_result.json().keys())
# dict_keys(['daily_quotes', 'pagination_key'])

pagination_key = daily_quotes_result.json()["pagination_key"]

# pagination_key も付与してもう一度リクエストする。
daily_quotes_url_2 = f"https://api.jquants.com/v1/prices/daily_quotes?date={date}&pagination_key={pagination_key}"
daily_quotes_result_2 = requests.get(daily_quotes_url_2, headers=headers)

# 今度は、pagination_keyは含まれていない。
print(daily_quotes_result_2.json().keys())
# dict_keys(['daily_quotes'])

# それぞれデータが得られている。
len(daily_quotes_result.json()["daily_quotes"]),  len(daily_quotes_result_2.json()["daily_quotes"])
# (4030, 312)

# それぞれ配列型のデータなので + で連結できる。
# DataFrame化までついでに行った。
df = pd.DataFrame(daily_quotes_result.json()["daily_quotes"]
                  + daily_quotes_result_2.json()["daily_quotes"])

print(len(df))
# 4342

1回目のリクエストでは、本当は4342件得られるはずだったデータのうち、4030件しか取得できてなかったことがわかりますね。そして、pagination_keyを合わせて送信することで、続きを取得できています。

上記のサンプルコードはわかりやすさ優先のため、2回で全部取得できると決め打ちしていますが、実際は2回目のリクエストでもpagination_keyが戻ってくる可能性があります。

そのため、実際の運用ではドキュメントのコードのようにpagination_keyがなくなるまでループするような実装にすると良いでしょう。

# 特定の日付の4本値を取得する
date = "2024-03-15"
daily_quotes_url = f"https://api.jquants.com/v1/prices/daily_quotes?date={date}"
headers = {"Authorization": f"Bearer {id_token}"}
daily_quotes_result = requests.get(daily_quotes_url, headers=headers)

# 1回目のレスポンスで得られたdata
data = daily_quotes_result.json()["daily_quotes"]

# pagination_keyが含まれている限りはループする。
while "pagination_key" in daily_quotes_result.json():
    pagination_key = daily_quotes_result.json()["pagination_key"]
    daily_quotes_url = f"https://api.jquants.com/v1/prices/daily_quotes?date={date}&pagination_key={pagination_key}"
    daily_quotes_result = requests.get(daily_quotes_url, headers=headers)
    # 得られたデータを連結する。
    data += daily_quotes_result.json()["daily_quotes"]


# データが揃っている。
print(len(data))
# 4342

これで、J-Quants APIのページング処理にも対応できました。

pandas.qcutでデータを分位数で離散化する

今回の記事ではpandasのqcutという関数を紹介します。
参考: pandas.qcut — pandas 2.2.1 ドキュメント

記事タイトルに書いていますが、これは分位数に基づいてデータを離散化する関数です。
実は以前、数値の区間で区切って離散化するpandas.cutというのを紹介したことがあるのですが、その仲間みたいなものですね。僕はつい最近までqcutを知りませんでしたが。
参考: pandasで数値データを区間ごとに区切って数える

cutでは数字の絶対値を基準に、0以上100未満、100以上200未満、みたいにデータを切り分けることができましたが、qcutでは分位数(パーセンタイル)を基準にデータを分けることができます。要するに、4つに分けるのであれば、25%以下、50%以下、75%以下、それより上、みたいにデータを区切り、各区切りには大体同じ件数のデータが分類されます。

cutだったら区間の幅が揃い、qcutだったら各区間に含まれるデータの件数が揃うというのが一番簡潔な説明ですね。

適当に乱数を使ってやってみましょう。ポアソン分布で200個ほどデータを作って、q_cutで5つのグループに分けてみます。

import pandas as pd
from scipy.stats import poisson  # テストデータ生成用


# λ=100のポアソン分布に従う乱数を200個生成
data = poisson(mu=100).rvs(size=200, random_state=0)
print(data[:10])  # 最初の10項表示
# [101 103  98  98 127 109 102  82  99  86]

# データと区切りたいグループの個数を指定して実行
out = pd.qcut(data, q=5)
# 各データがそれが含まれる区間お
print(out)
"""
[(97.0, 102.0], (102.0, 107.0], (97.0, 102.0], (97.0, 102.0], (107.0, 127.0], ..., (102.0, 107.0], (102.0, 107.0], (92.0, 97.0], (92.0, 97.0], (78.999, 92.0]]
Length: 200
Categories (5, interval[float64, right]): [(78.999, 92.0] < (92.0, 97.0] < (97.0, 102.0] < (102.0, 107.0] < (107.0, 127.0]]
"""

データの先頭の方と、あと、結果をprintして表示されたやつを上のコードに出しました。Categories として5つの区間が表示されていますが、「それぞれのデータがどの区間に含まれているのか」に変換されたものが得られていますね。例えば最初のデータは101ですが、これは区間(97, 102] に含まれます。

区間にラベルをつけることもできます。低い方からL1, L2, L3 みたいにつけていく場合はlabel引数にqで指定した数と同じ要素数の配列を渡して実現します。(今回文字列でサンプル作っていますが、数値をラベルにすることもできます。)

# ラベルを指定する
out = pd.qcut(data, q=5, labels=["L1", "L2", "L3", "L4", "L5"])
print(out)
"""
['L3', 'L4', 'L3', 'L3', 'L5', ..., 'L4', 'L4', 'L2', 'L2', 'L1']
Length: 200
Categories (5, object): ['L1' < 'L2' < 'L3' < 'L4' < 'L5']
"""

変換後のデータとして扱いやすそうな形で結果が得られました。

ただ、それぞれのラベルの区間がこれだとわからないですね。区間の情報を別途得る必要があるのでその場合はretbins引数にTrueを渡して、結果を受け取るときにもう一個変数を用意して受け取ることで、区切り位置の譲歩を得ることもできます。もちろん、labelsは使わずに、retbinsだけ指定することもできますよ。

# ラベルを指定する
out, bins = pd.qcut(data, q=5, labels=["L1", "L2", "L3", "L4", "L5"], retbins=True)
print(out)
"""
['L3', 'L4', 'L3', 'L3', 'L5', ..., 'L4', 'L4', 'L2', 'L2', 'L1']
Length: 200
Categories (5, object): ['L1' < 'L2' < 'L3' < 'L4' < 'L5']
"""
# 区切り位置の情報
print(bins)
# [ 79.  92.  97. 102. 107. 127.]

最後に注意です。qcutを使うと連続値のデータは大体同じ個数ずつに分けてくれることが多くそれが目的で使うことが多くなるのですが、今回の例のように整数値など離散な値しか取らない場合はそうでもなくなってきます。今回乱数で発生したデータはちょうど区切り位置の107が10個も混ざってた等々の事情で、ちょっとだけ偏りが出ています。実際に使う場合はこのあたりの結果もよく注意してみてください。

print(out.value_counts())
"""
L1    44
L2    39
L3    39
L4    41
L5    37
dtype: int64
"""

np.vectorizeで関数をベクトル化する

NumPyやScyPyの関数って非常に便利で、NumPy配列(要するにArray)を渡すと空気を読んでその渡したデータの各要素に関数を適用してNumPy配列で結果を返してくれたりします。

自分で定義した関数でもNumPyやSciPyの関数の組み合わせで作った関数であれば結構そのように動いてくれるのですが、文字列操作が入ったりif文による分岐等があると必ずしもそうはならず、スカラー値を受け取ってスカラー値を返すだけの関数になることがあります。

そのような関数を、手軽にベクトルか対応することができる方法があるのでこの記事で紹介します。

それが、記事タイトルのnp.vectorizeです。

ドキュメント: numpy.vectorize — NumPy v1.26 Manual

関数を渡すと戻り値で新しい関数オブジェクトが帰ってきてそれがベクトル対応(配列対応)しています。

基本的な使い方

数学関数だと特にArrayを渡すと元々期待通り動いたりするので、少々無理矢理な例ですが文字列操作の関数を作ってお見せします。これは数値を1個受けとって、その数値に、「回目」っていいう単位をつけて返すだけの関数です。普通に実験、そのまま配列渡してみる、ベクとライズして配列を渡してみる、の3パターンやってみました。

import numpy as np


# 数値に単位をつける関数を実装
def number_format(n):
    return f"{n}回目"


# 数値を渡すと想定通り動く
print(number_format(5))
# 5回目

# 配列を渡すと配列を一個の値とみなして文字列化して単位をつけてしまう。
print(number_format([1, 2, 3]))
# [1, 2, 3]回目

# ベクトル化した関数を作る
number_format_vec = np.vectorize(number_format)

# それに配列を渡すと配列の各要素に元の関数を適用してくれる。
print(number_format_vec([1, 2, 3]))
# ['1回目' '2回目' '3回目']

# Array型もタプルもいける
print(number_format_vec(np.array([1, 2, 3])))
# ['1回目' '2回目' '3回目']
print(number_format_vec((1, 2, 3)))
# ['1回目' '2回目' '3回目']

# もちろん、内包表記で同じことをすることは可能(ただし、この結果はlist)
print([number_format(n) for n in [1, 2, 3]])
# ['1回目', '2回目', '3回目']

ベクトル化した関数を1回しか使わないなら内包表記で済ましちゃっていいんじゃないかな、と思うのですが、何度も利用したい関数であればnp.vectorizeを使うと言う選択肢もあるのかな、と思います。

注意点

NumPyやSciPyで実装されている関数群って並列処理できる部分は並列処理するような賢い実装になっていることがありますが、この np.vectorize はそこまで気が利いたものではありません。どうやら単純にfor文で順次処理するようになるだけらしいので処理の高速化等の効果はありません。ドキュメントにも利便性のためのもので、パフォーマンスのため使うようなものではなく、for loop回してるだけだって書いてありますね。

そのため、本当に頻繁に大規模なベクトルを処理する関数なのであれば別の方法で対応させる必要があるでしょう。

もう一点、細かいですが戻り値がNumPyのarrayであることも注意が必要ですね。と言ってもこれは便利に感じることが多いですが。内包表記であればlistで結果が得られますがvectorizeするとlist渡してもlistではなくarrayで帰ってきます。

引数を複数受け取る関数の場合

この np.vectorize は引数を複数受け取る関数にも対応しています。ドキュメントのサンプルもa, b の2変数受け取っていますしね。一応その例も見ておきましょう。年と月の数値を受け取って何年何月、という文字列返す関数でやってみます。

def month_str(year, month):
    return (f"{year}年{month}月")


month_str_vec = np.vectorize(month_str)

# 元の関数はyear, monthは1個ずつしか値を受け取れない
print(month_str([2020, 2023, 2026], [1, 4, 7]))
# [2020, 2023, 2026]年[1, 4, 7]月

# ベクトル化すると複数ペアをまとめて処理できる。
print(month_str_vec([2020, 2023, 2026], [1, 4, 7]))
# ['2020年1月' '2023年4月' '2026年7月']

# 片方は配列で、片方はスカラーというパターンにも対応する
print(month_str_vec([2020, 2023, 2026], 1))
# ['2020年1月' '2023年1月' '2026年1月']

さいごに

以上が手軽に関数をベクトル化する方法でした。まぁ、内包表記もあればmapを使うやり方もあるのでこれが必須というわけではないのですがいい感じに動く関数を手軽に作る方法として頭の片隅に置いておくと使う場面はあるんじゃないかなと思います。

ちなみに、関数を定義した直後にベクトル化した関数で元の関数名を上書きしておくと、最初っからベクトル化した関数を宣言したのと同じように使えますよ。

def func(x):
    # 何かの処理


func = np.vectorize(func)
# 以降に呼び出されるfuncはベクトル対応した関数。

Macで特定のポートが使用中かどうか調べる

Jupyter LabをローカルのMacで使っており、それを起動するバッチを自作して使っているのですが、たまにすでに起動してるのにそのバッチを実行してしまってJupyterのプロセスを2重に立ち上げてしまうことが度々ありました。

2個立ち上げっぱなしとなると気持ち悪いのでPIDを調べてkillする必要がありやらかすとちょっと面倒なミスです。

最初、この対策としてプロセスをgrepしてjupyterがすでに立ち上がってたらバッチを中断して起動コマンドを叩かない、みないな処理を入れて対応しようとも試みていたのですが、grep jupyter 自身がそのgrepにヒットするというなかなか困った挙動をしうまくいかずに放置していました。

しかしその後、プロセス一覧ではなくポートが使用中かどうかで判断すれば良いと思いついたので実装してみました。

利用するのは、lsof というコマンドです。これはポートに限らず、システムで開いているファイルの一覧を表示するコマンドです。名前も List Open Files の略だそうです。

今回はポートの情報だけ知れたらいいので、 -i 引数をつけます。そして、さらにコロンを打ってポート番号を指定すると、そのポートが使用されていれば使用しているプロセスの一覧が各種情報とともにずらずらと表示されます。

Jupyter Labは デフォルトでは 8888番ポートで動作するので、次のようなコマンドで確認できます。

$ lsof -i :8888
# 8888番ポートが使われていなければ何も表示されない。

# Jupyter Labが起動していると各プロセスが結果に出る。
COMMAND    PID   USER   FD   TYPE            DEVICE SIZE/OFF NODE NAME
Google     808 yutaro   21u  IPv6 0x12572c541560807      0t0  TCP localhost:49319->localhost:ddi-tcp-1 (ESTABLISHED)
Google     808 yutaro   34u  IPv4 0x12572d3a3516cff      0t0  TCP localhost:49234->localhost:ddi-tcp-1 (ESTABLISHED)
Google     808 yutaro   37u  IPv6 0x12572c541540007      0t0  TCP localhost:49253->localhost:ddi-tcp-1 (ESTABLISHED)
Google     808 yutaro   80u  IPv6 0x12572c54155b807      0t0  TCP localhost:49213->localhost:ddi-tcp-1 (ESTABLISHED)
Google     808 yutaro  108u  IPv6 0x12572c541546007      0t0  TCP localhost:49216->localhost:ddi-tcp-1 (ESTABLISHED)
python3.1 1096 yutaro   11u  IPv4 0x12572d3a351ba6f      0t0  TCP localhost:ddi-tcp-1 (LISTEN)
python3.1 1096 yutaro   12u  IPv6 0x12572c54154a007      0t0  TCP localhost:ddi-tcp-1 (LISTEN)
python3.1 1096 yutaro   13u  IPv6 0x12572c54155e807      0t0  TCP localhost:ddi-tcp-1->localhost:49319 (ESTABLISHED)
python3.1 1096 yutaro   14u  IPv6 0x12572c54155c007      0t0  TCP localhost:ddi-tcp-1->localhost:49213 (ESTABLISHED)
python3.1 1096 yutaro   16u  IPv6 0x12572c541548007      0t0  TCP localhost:ddi-tcp-1->localhost:49216 (ESTABLISHED)
python3.1 1096 yutaro   37u  IPv4 0x12572d3a351780f      0t0  TCP localhost:ddi-tcp-1->localhost:49234 (ESTABLISHED)
python3.1 1096 yutaro   50u  IPv6 0x12572c54153e807      0t0  TCP localhost:ddi-tcp-1->localhost:49253 (ESTABLISHED)

Python は3.11を使っているのに、COMMAND名が途中で切り捨てられて3.1って表示されていますが、まぁ、Pythonで何か動いているのはわかりますね。

あとは、Jupyter Labの起動バッチで、このコマンドを動かして処理を分岐させればOKです。

(自分はnohupコマンドの結果を専用のログファイルにログを書き出しているのですが、不要なら/dev/null あたりに捨てても良いと思います。)

if 文の条件文の中で出力を /dev/nullに捨てているのは出力をすっきりさせるためです。捨てない場合は普通に画面にlsofコマンドの結果も表示されます。

#!/usr/bin/env zsh
if lsof -i:8888 > /dev/null; then
    echo "8888番ポートは既に使用されています。"
    exit 1
else
    cd ~
    nohup jupyter lab >> ~/Documents/log/jupyter_lab.log 2>&1 &!
    cd -
fi

cd ~ はホームディレクトリに移動する、 cd – は元のディレクトリに戻る、です。
どこで打ってもホームディレクトリでJupyter Labが起動されるように、そしてそのあとは元のディレクトリに帰るようにしています。

この記事を書いた後に気づいきましたが、このlsofコマンドは以前こちらの記事でも使っていましたね。
参考: sshtunnel を使って踏み台サーバー経由でDB接続

コマンドでMacbookのスリープを一時的に抑制する

バッテリーの持ちとセキュリティ関連の理由により、一定時間アイドル状態になったらスリープする設定でMacを使っている方は多いと思います。
僕もそうです。

ただ、機械学習のモデルを学習している場合や、ベイズモデルでMCMCのサンプリングをやっている時など、しばらく作業の待ち時間が発生することがあり、その脇で本を世だりして時間を潰しているといつの間にかPCがスリープしてしまうということがあります。

普通にシステム環境設定からスリープまでの時間の設定を変えて、作業後に戻せばいいだけの話なのですが、それはやや手間です。

しかし、最近コマンドで一時的にスリープを抑制できることを知ったのでそれを紹介します。

コマンド名は caffeinate です。名前はコーヒーのカフェインが由来だとか。面白いですね。

細かいオプションが複数あり、それらを組み合わせて使います。

-d ・・・ ディスプレイのスリープを抑制
-i ・・・ システムのスリープを抑制
-m ・・・ ディスクのスリープを抑制
-s ・・・ AC電源で動作している場合にシステムのスリープを抑制
-u ・・・ -tとセットで利用し、ユーザーがアクティブであることを宣言する
-t ・・・ -uの時間を秒数で指定する

例えば30分間くらいスリープしたくないのであればこんな感じですね。

$ caffeinate -d -u -t 1800

-t で時間を指定するとその時間経過時に勝手に終了しますが、そうでない場合は Ctrl + d でコマンドを中断するまで抑制してくれます。

そもそも一定時間でスリープするような設定にしている事情(おそらくセキュリティ的な話)があると思うので、個人的には-tで時間を指定して使うのがお勧めです。

scipy.statsで確率分布を自作する

SciPyのstatsモジュールには非常に多くの確率分布が定義されています。
参考: Statistical functions (scipy.stats) — SciPy v1.12.0 Manual

ほとんどの用途はこれらを利用すると事足りるのですが、自分で定義した確率分布を扱うこともあり、scipyのstatsに用意されているような確率分布と同じようにcdfとかのメソッドを使いたいなと思うことがあります。

そのような場合、SciPyではscipy.stats.rv_continuousを継承してpdfメソッドを定義することで確率分布を定義でき、そうすれば(計算量や時間の問題などありますが)cdf等々のメソッドが使えるようになります。

参考: scipy.stats.rv_continuous — SciPy v1.12.0 Manual

実はSciPyの元々用意されている確率分布たちもこのrv_continuousを継承して作られているので、それらのソースコードを読むと使い方の参考になります。既存の確率分布等はpdfだけでなく、cdf、sf、ppf、期待値や分散などの統計量などが明示的に実装されていて計算効率よく使えるようになっています。

しかし先ほども述べた通り、pdf以外は必須ではなく実装されていなければpdfから数値的に計算してくれます。(ただ、予想外のところで計算が終わらなかったり精度が不十分だったりといった問題が起きるので可能な限り実装することをお勧めします。)

確率分布クラスの自作方法

さて、それでは実際にやってみましょう。既存に存在しない例が良いと思うので、正規分布を2個組み合わせた、山が2個ある分布を作ってみます。

クラスを継承して _pdf (pdfではない) メソッドをオーバーライドする形で実装します。

import numpy as np
from scipy.stats import rv_continuous
import matplotlib.pyplot as plt


# 確率密度関数の定義
def my_pdf(x):
    curve1 = np.exp(-(x+3)**2 / 2) / np.sqrt(2 * np.pi)
    curve2 = np.exp(-(x-3)**2 / 2) / np.sqrt(2 * np.pi)

    return (curve1 + curve2)/2


# rv_continuousクラスのサブクラス化
class MyDistribution(rv_continuous):
    def _pdf(self, x):
        return my_pdf(x)


# 分布のインスタンス化
my_dist = MyDistribution(name='my_dist')


# 可視化
x = np.linspace(-6, 6, 121)
fig = plt.figure(facecolor="w", figsize=(8, 8))
ax = fig.add_subplot(2, 2, 1, title="pdf")
ax.plot(x, my_dist.pdf(x))

ax = fig.add_subplot(2, 2, 2, title="cdf")
ax.plot(x, my_dist.cdf(x))

ax = fig.add_subplot(2, 2, 3, title="sf")
ax.plot(x, my_dist.sf(x))

ax = fig.add_subplot(2, 2, 4, title="logpdf")
ax.plot(x, my_dist.logpdf(x))

plt.show()

結果がこちらです。

実装したのはpdfだけでしたが、cdfやsf、logpdfなども動いていますね。
繰り返しになりますが、これらは数値計算されたものなので十分注意して扱ってください。特に精度と計算量はもちろんですが、指数関数等もからむので、極端な値を入れたら計算途中でオーバーフローが起きたりもします。重要なものは_pdfと同様に_cdfなどとして明示的に定義した方が良いでしょう。

たとえば、先ほどの分布は分散の計算に少し時間がかかります。

%%time
my_dist.var()
"""
CPU times: user 8.02 s, sys: 146 ms, total: 8.16 s
Wall time: 8.06 s
10.000000000076088
"""

確率密度関数を用意するときの注意点

確率密度関数を自分で作りましたが、この内容に対してのバリデーションなどは用意されていないので自分で責任を持って管理する必要があります。

たとえば、変数xの定義域で正であること、定義域全体で積分したら1になることなどは事前に確認しておく必要があります。(確率密度関数の定義を満たしてないととcdfなど他のメソッドが変な動きをします)

確率分布の台が有限の場合

もう一点、先ほどの2山の分布はxが$-\infty$から$\infty$の値を取る分布でしたが、確率分布の中にはそうではないものもあります。正の範囲でしか定義されない対数正規分布や指数分布、有限区間でしか定義されないベータ分布などですね。

これらについては、「分布のインスタンス化」を行うタイミングで台の下限上限をそれぞれa, b という引数で行います。省略した方は$-\infty$もしくは$\infty$となります。

確率密度関数が上に凸な放物線で定義された分布でやるとこんな感じです。

def my_pdf2(x):
    return 6*x*(1-x)


class MyDistribution2(rv_continuous):
    def _pdf(self, x):
        return my_pdf2(x)


# 分布のインスタンス化
my_dist2 = MyDistribution2(a=0, b=1, name='my_dist2')

以上が、SciPyで連続確率分布インスタンスを自作する方法でした。

SciPyで重積分

もう結構古い記事なのですが、以前SciPyで定積分をやる方法を記事にしたことがあります。
参考: scipyで定積分

最近、2変数関数の積分をやる機会があったのでこの機会に重積分をSciPyで行う方法を紹介します。SciPyのintegrateモジュールには、重積分用の関数が複数あります。
dblquad (2変数関数の定積分)
tplquad (3変数関数の定積分)
nquad (一般のn変数の定積分)

dblquadの使い方

順番に説明していきます。まずは2重積分のdblquadです。関数の定義は次のようになっています。
scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

必須なのは、積分対象のfunc, 外側の積分区間のa, b、そして内側の積分区間を示す、gfun, hfunです。

gfunとhfunは名前からわかる通り、定数ではなく関数です。これにより内側の積分の積分区間を変数にすることができます。つまり以下のような積分区間の積分ができます。
$$\int_0^1\int_0^y xy \,dxdy$$

例に挙げたのでこれを実装してみましょう。ちなみに解は$1/8=0.125$です。funcの定義は、内側の関数を第1引数にする必要があるので注意してください。

from scipy.integrate import dblquad


def f1(x, y):
    return x*y


def x0(y):
    return 0


def x1(y):
    return y


print(dblquad(f1, 0, 1, x0, x1))
# (0.125, 5.515032205777789e-15)

想定通りですね。積分結果と推定誤差が返ってくるのは1変数の積分と同様です。

内側の積分区間も定数から定数までだよ、要するに長方形領域で積分したいよ、って場合はgfun, hfun に定数を返す関数を返してください。

tplquadの使い方

続いて、3変数向けのtplquadです。これもdblquadとかなり似てる感じで使えます。積分変数が一個増えているので上限加減の指定がもう一個ある感じです。
scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

たとえば次の積分をやってみましょう。

$$\int_0^1\int_0^z\int_0^{y+z} xyz \,dxdydz.$$

ちなみに答えは$17/144=0.11805555…$となるはずです。

引数の順番に注意が必要なので慎重にコーディングしてください。

from scipy.integrate import tplquad


def f2(x, y, z):
    return x*y*z


def x0(y, z):
    return 0


def x1(y, z):
    return y + z


def y0(z):
    return 0


def y1(z):
    return z


print(tplquad(f2, 0, 1, y0, y1,  x0, x1))
# (0.11805555555555557, 2.1916761217856673e-14)

バッチリですね。

nquadの使い方

最後に一般のn変数を積分できるnquadの使い方を紹介します。

引数の形式が先ほどの二つと少し違います。
scipy.integrate.nquad(func, ranges, args=None, opts=None, full_output=False)

funcにn変数関数を渡して、rangesに積分区間を渡すことになります。rangesは配列で、1変数目から順番に区間の下限上限の2値の配列を格納しておけば良いです。また、ここにも一応関数を使うことはできます。

これはシンプルな例で、定数関数1を超立方体区間で積分してみました。

from scipy.integrate import nquad


def f3(w, x, y, z):
    return 1

print(nquad(f3, [[-2, 2], [-2, 2], [-2, 2], [-2, 2]]))
# (256.0, 2.8421709430404007e-12)

$4^4=256$になりましたね。

ここで急にシンプルな例を出したのには事情がありまして、変数の数が多くなるとやはり積分は困難なようで、ちょっと複雑な例になると規定の反復回数をこなしても必要な精度に届かずWarningが出たりするケースが多々あります。

どうしても計算したい場合は limit パラメーター等をいじっての対応になりますのでドキュメントを参照しながら調整してみてください。
(僕も実運用で必要になったら改めて調査して紹介しようと思います。)

Pythonの関数から一部の引数を固定して新しい関数を作る

Pythonの多くのライブラリの様々な関数が非常に汎用的に使えるように作られているので多くの引数を受け取れるようになっています。しかし、そのほとんどの引数を固定して1変数関数として使いたいなぁと思うようなことがあります。PandasのDataFrameのapplyなど関数を引数として受け取る関数に渡す場合等ですね。
また、大量にある引数のほとんどを固定して一部だけ変えながら何度も実行する、といった場面も考えられます。

lambda式などを作ってラップした新しい関数を実装してもいいのですが、 functoolsという標準ライブラリにその専用のpartial というメソッドが用意されています。
参考: functools.partial(func/*args**keywords)

このpartialを使うと、引数の一部を固定した引数の少ない新しい関数を作ってくれます。

一個目の引数に元になる関数を渡し、2個目以降の引数に渡したものが、元の関数の固定引数として使われます。keyword引数で渡せばそのkeyword引数が固定されます。

一引数の固定の方は先頭から順番に固定されるので注意してください。つまり2番目以降の引数を固定したい場合はそれらはキーワード引数として指定する必要があります。

サンプル

引数を順番に表示するだけの単純な関数を作ってやってみましょう。

from functools import partial


# 3つの引数を表示するだけの関数
def sample_func(a, b, c):
    print("a=", a)
    print("b=", b)
    print("c=", c)


# テスト実行
sample_func(1, 2, 3)
"""
a= 1
b= 2
c= 3
"""

# a = 10, b = 20 を固定した新しい関数が作られる。
partial_f = partial(sample_func, 10, 20)


# 3個目の引数 c = 50だけ渡して実行できる。
partial_f(50)
"""
a= 10
b= 20
c= 50
"""

# キーワード引数で固定することもできる。
partial_f2 = partial(sample_func, a=100, c=200)

# b の値だけ渡して実行できる
partial_f2(b=-5)
"""
a= 100
b= -5
c= 200
"""

キーワード引数を固定した関数を、位置引数で使う場合は注意が必要です。
たとえば、次のようにaを固定して生成した関数に、残り2個の引数を位置引数で渡すと、aを2回渡した扱いになってエラーが起きます。

# aを固定
partial_f3 = partial(sample_func, a=1)

# bとcのつもりで残り2個の引数を渡すとエラー
try:
    partial_f3(2, 3)
except Exception as e:
    print(e)
# sample_func() got multiple values for argument 'a'

# bとcもキーワード引数で渡す。
partial_f3(b=2, c=3)
"""
a= 1
b= 2
c= 3
"""

まとめ

ほぼ小ネタのような内容でしたが、自作関数をベースに一部の振る舞いを固定した簡易的な関数を作るとか、apply等の1変数関数を受け取るメソッドに渡したいとかそういう場面で役に立つことがあるテクニックとしてpartialを紹介しました。

scipyのstats配下の各種メソッドであれば、それぞれがパラメーターを固定するfrozenメソッドを持ってるとか、引数が多いなら引数を辞書にまとめて**(アスタリスク2個)で展開すればいいとか、ラップした関数を自分で実装したらいいとか、代用手段も多いのですが、partialを使うとその辺の記述がシンプルになるので機会があれば使ってみてください。

SciPyでニュートン法を利用する

前回の記事の二分法に続いて、もう一つ求根アルゴリズムを紹介します。
参考: 二分法を用いて関数の根を求める

今回紹介するのはニュートン法です。これは微分可能な関数$f(x)$の根を求めることができるアルゴリズムです。
参考: ニュートン法 – Wikipedia

詳しい説明は上記のWikipediaにあるので、ざっくりと概要を説明します。

この方法の背景にあるのは、滑らかな関数をある点の近くだけ着目してみるとほぼ直線になり、接線で近似できるということをベースのアイデアにしています。

つまり、微分可能な関数$f$があって$f(x)=0$だとします。その根$x$の近くに点$x_0$を取ると、$x_0$の近くでは、$f$と$f$の接線ってかなり近いよね、それなら$f$の根と$f$の$x_0$における接線の根って近いよね、っていうのが基本的なアイデアです。

関数$f$の$f(x_0)$における接戦は次の式で書けます。

$$y=f'(x_0)(x-x_0)+f(x_0).$$

$f(x)=0$は解けない場合でも、この接線の根は容易に算出することができ、

$$x_0-\frac{f(x_0)}{f'(x_0)}$$

と求まります。

この値は元の$x_0$よりも真の根に近いことが期待され、これをもう一回$x_0$とおいて同じ操作を繰り返せば真の根にたどり着く、というのがニュートン法です。

Wikipediaから画像拝借しますが、図で見るとイメージしやすいですね。

注意しないといけないのは、初期値$x_0$は真の根$x$の十分近くに取らないといけない点です。十分近くを見れば関数をその接線で近侍できるよね、というのがアイデアの前提なので、根が近くになかったらその前提が崩れてしまいこのアルゴリズムは真の根に収束しなくなってしまいます。

ニュートン法のメリットとデメリット

先に紹介した二分法と比べて、ニュートン法のメリットデメリットを説明していきます。

1番のメリットは収束の速さです。二分法に比べてより少ない計算回数で効率的に会を探索することができます。

また、初期値として与える点が1点だけで良いというのもメリットです。二分法の場合は初期値は区間で設定する必要がありましたからね。

その一方で複数のデメリットもあります。実装していて一番不便に感じるのはその関数だけでなく微分も必要ということでしょうか。もちろん微分不可能な関数ではニュートン法は使えません。

また、初期値が真の解に十分近くない場合や、微分した値が$0$に近い場合、うまく収束せずにアルゴリズムが失敗してしまう、という点も大きなデメリットです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。newtonという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’newton’)を指定して使うことになります。二分法と同じですね。

参考:
scipy.optimize.newton — SciPy v1.12.0 Manual
root_scalar(method=’newton’) — SciPy v1.12.0 Manual

二分法の時と同じように、$\sin$関数の根$\pi$を探索させてみましょう。微分は$\cos$なのでこれを使います。

from scipy import optimize
import numpy as np


root1 = optimize.newton(np.sin, x0=3, fprime=np.cos)
print(root1)
# 3.141592653589793

root_result = optimize.root_scalar(np.sin, method="newton", x0=3, fprime=np.cos)
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 6
     iterations: 3
           root: 3.141592653589793
"""

print(root_result.root)
# 3.141592653589793

簡単ですね。

注目するのは、iterationsの部分です。たった3回のイテレーションで収束していて、関数が実行されたのは、fとfの微分合わせて6回だけです。
二分法の時は39回もイテレーションが必要だったのと大違いです。そして実はこの例では解の精度もニュートン法の方が高くなっています。

ニュートン法が失敗する例

初期値が真の解の近くにないと失敗するという話がありましたのでそちらも見ておきます。

例えば、タンジェントの逆関数、$\arctan$で試してみましょう。(sin, cosは根が無限にあって、根から遠い実数を用意できないので関数を変えます。)

$\arctan(x)$の微分は$\frac{1}{1+x^2}$です。

やってみました。

def f(x):
    return np.arctan(x)


def fprime(x):
    return 1/(1+x**2)


# 初期値が1なら収束する
root_result_1 = optimize.root_scalar(np.sin, method="newton", x0=1, fprime=fprime)
print(root_result_1)
"""
      converged: True
           flag: 'converged'
 function_calls: 12
     iterations: 6
           root: 0.0
"""

# 初期値が2だと失敗し、結果のflagが'convergence error'になる。
root_result_2 = optimize.root_scalar(np.sin, method="newton", x0=2, fprime=fprime)
print(root_result_2)
"""
      converged: False
           flag: 'convergence error'
 function_calls: 100
     iterations: 50
           root: 1.854706857103781
"""


# optimize.newton の方だと例外が上がる。
try:
    optimize.newton(f, fprime=fprime, x0=2)
except Exception as e:
    print(e)
# Derivative was zero. Failed to converge after 10 iterations, value is -6.999943395317963e+168.

失敗した時の振る舞いがそれぞれ違うので、どちらのコードを使うかで注意深く扱う必要がありますね。optimize.root_scalarはコード自体は正常に終了しますがフラグが立ち、optimize.newtonの方は例外があがります。

ちなみに、エラーの中で出てくるvalue の値、 -6.999943395317963e+168 は次のように自分でニュートン法を実装しても同じ値が出て来ます。

x0 = 2  # 初期値
for i in range(12):
    x0 = x0 - f(x0)/fprime(x0)
    print(i+1, "回目: x0=", x0)

"""
1 回目: x0= -3.535743588970453
2 回目: x0= 13.950959086927496
3 回目: x0= -279.34406653361754
4 回目: x0= 122016.9989179547
5 回目: x0= -23386004197.933937
6 回目: x0= 8.590766671950415e+20
7 回目: x0= -1.1592676698907411e+42
8 回目: x0= 2.110995587611039e+84
9 回目: x0= -6.999943395317963e+168
10 回目: x0= inf
11 回目: x0= nan
12 回目: x0= nan
"""

絶対値が大きくなり続けていて全く収束に向かっていないのがわかりますね。

まとめ

元の関数だけではなく導関数も必要だったり、初期値の設定段階である程度解の目星をつけておかないといけないなどのデメリットはありますが、速度や精度の面で優秀でしかもロジックもわかりやすい手法なので、何か機会があればニュートン法の活用を検討してみてください。