静かなる名辞

pythonとプログラミングのこと

2019/03/22:TechAcademyがteratailの質問・回答を盗用していた件
2019/03/26:TechAcademy盗用事件 公式発表と深まる疑念



【python】print関数を使いこなそう

ぼくたちは本当のprintを知らない

 pythonのprint関数については、たかがprintと思っている人も多いと思いますが、しかしオプションをすべて言える人はあまりいないと思います。把握しておくと出力の細かい制御をしたいとき役立ちます。

 そこで、printの使いこなしについて書きます。なお、念のために先に断っておくと、新し目のpython3のprint関数の話です。python2のprint文とはまったく異なります。

printのオプション

 ドキュメントによると以下の引数を取ります。

print(*objects, sep=' ', end='\n', file=sys.stdout, flush=False)

組み込み関数 — Python 3.7.3rc1 ドキュメント

 *objectsは任意個のオブジェクトをprintできることを示します。sep=' 'は、オブジェクトの区切り文字が半角スペースになるということです。これは*objectsに渡したオブジェクトの数-1個表示されます。また、end='\n'は行末に何を置くかを指定しており、デフォルトは改行です。これは行末に1つだけ置かれます。

 特殊なのはfile引数とflush引数です。file引数は出力ストリームに対応し、デフォルトのsys.stdoutはいわゆる標準出力です。任意のファイルオブジェクトを与えられるので、ファイルに書き込むといった用途に使えます。flushはいわゆる出力バッファのフラッシュの挙動になります。何もしなければ標準出力は改行時にフラッシュされますが、これをTrueにするとprintを呼ぶたびにフラッシュされるようになります。

sepの使い方

 何も指定しないと半角スペースになるので、次のような見慣れた挙動になります。

>>> print(1, 2, 3)
1 2 3

 空の文字列にすることで間を詰めることが可能です。

>>> print(1, 2, 3, sep="")
123

 あるいは好きな文字列を渡すことも可能です。

>>> print(1, 2, 3, sep=",")
1,2,3
>>> print(1, 2, 3, sep="***")
1***2***3

 リストの要素を適当な区切り文字で区切って表示したいというような場合、str.joinで頑張るより、リストのアンパックとprintのsep引数を組み合わせて使った方がスマートだったりします。

>>> l = [1, 2, 3]
>>> print(",".join(map(str, l)))  # こちらの場合は型変換も必要
1,2,3
>>> print(*l, sep=",")  # これだけ
1,2,3

 裏を返せば表示形式はすべて__str__任せになるので、細かい出力フォーマットは指定できないということでもあるのですが、簡易的に使いたい場合はsepの方が便利です。

endの使い方

 改行させたくないときにendを空文字列にする、というのがよくあるユースケースです。

>>> for x in range(4):
...     print(x)
... 
0
1
2
3
>>> for x in range(4):
...     print(x, end="")
... 
0123

 最後に一回改行するために、もう一回printを呼ぶと良いと思います。

 プログレスバー的な進捗表示を簡単にやりたいとき重宝しますが、それっぽく見せるためには後述のflush引数も併用する必要があります。

fileの使い方

 fileはたとえばファイル出力に使えます。

>>> with open("hoge.txt", "w") as f:
...     for x in range(4):
...         print(x, file=f)
... 
>>> (ctrl-dでpythonを終了)
$ cat hoge.txt
0
1
2
3

 その気になればcsvファイルの文字列を手動で作って吐き出すといった用途に使えます。あるいはソケットを叩いたりするのにも使えると思います。ただ、濫用しすぎないことをおすすめします(通常は他に良い方法がある)。

flushの使い方

 flushは前述した通りendと組み合わせると便利です。

>>> import time
>>> for _ in range(10):
...     print("*", end="")
...     time.sleep(0.5)
... 
**********

 進捗表示をやりたくてこんなコードを書く人も多いと思いますが、これは期待通り動作しないと思います。確かに最終的な結果は同じでも、5秒ほどしてから出力が一気に出てくるという挙動になるはずです。これはsys.stdoutが改行されるまでフラッシュされないからです。

 flush=Trueで解決します。

>>> import time
>>> for _ in range(10):
...     print("*", end="", flush=True)
...     time.sleep(0.5)
... 
**********

 これでおよそ0.5秒おきに1つずつ表示されるはずです。

本当のprintを理解したあとの感想

 print一つとっても、意外といろいろな使い方があるということがわかるかと思います。出力の制御にけっこう使えます。

 こういう仕様はすべて公式ドキュメントに書いてあります。ドキュメントに普段から慣れ親しんでおくと、スムーズにコーディングできる、行き詰まったとき打開策が思い浮かぶ、ライバルと差が付けられる、といった利益があります。ドキュメントを読むのに慣れていない方は、printなど単純なものから読んでみると楽しいと思います。

【python】sklearnのRandomizedSearchCVを使ってみる

はじめに

 RandomizedSearchCVなるものがあるということを知ったので、使ってみます。うまく使うとグリッドサーチよりよい結果を生むかもしれないということです。

sklearn.model_selection.RandomizedSearchCV — scikit-learn 0.20.3 documentation

比較実験

 とりあえず、先に使い慣れたグリッドサーチでやってみます。digitsデータをSVMで分類するという、100回くらい見た気がするネタです。

コード

import time

import pandas as pd

from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
from sklearn.datasets import load_digits

def main():
    digits = load_digits()
    svm = SVC()
    params = {"C":[0.1, 1, 10], "gamma":[0.001, 0.01, 0.1]}
    clf = GridSearchCV(svm, params, cv=5, iid=True,
                       return_train_score=False)
    t1 = time.time()
    clf.fit(digits.data, digits.target)
    t2 = time.time()
    print("{:.2f} 秒かかった".format(t2 - t1))
    result_df = pd.DataFrame(clf.cv_results_)
    result_df.sort_values(
        by="rank_test_score", inplace=True)
    print(result_df[["rank_test_score", 
                     "params", 
                     "mean_test_score"]])

if __name__ == "__main__":
    main()

結果

16.67 秒かかった
   rank_test_score                      params  mean_test_score
6                1   {'C': 10, 'gamma': 0.001}         0.972732
3                2    {'C': 1, 'gamma': 0.001}         0.971619
0                3  {'C': 0.1, 'gamma': 0.001}         0.943239
7                4    {'C': 10, 'gamma': 0.01}         0.705620
4                5     {'C': 1, 'gamma': 0.01}         0.695047
1                6   {'C': 0.1, 'gamma': 0.01}         0.111297
5                7      {'C': 1, 'gamma': 0.1}         0.104062
8                7     {'C': 10, 'gamma': 0.1}         0.104062
2                9    {'C': 0.1, 'gamma': 0.1}         0.102393

 まあ、これはこんなもんでしょう。しっかり最適パラメータを探せている気がします。3*3=9通りのグリッドサーチで、計算には16.67秒かかったとのことです。

RandomizedSearchCVにしてみる

 RandomizedSearchCVの特色は、scipyで作れる確率分布のオブジェクトを渡せることです。パラメータのリストを渡すことも可能ですが、それだと特色を活かした使い方にはなりません。

 scipyで確率分布のオブジェクトを作る方法については、以前の記事で説明したのでこちらを見てください。

www.haya-programming.com

 ここで期待されている「確率分布のオブジェクト」はrvs(分布に従うサンプルを得るメソッド)さえ使えれば何でも良いです。連続分布も離散分布も渡せます。なんなら自作でも行けると思います。

 といってみても、パラメータ設定を確率分布に落とし込むまでは少し頭の体操が要ると思います。

 次のように考えます。たとえば、SVMでパラメータチューニングするなら、とりあえずこんな感じでざっくり見る人が多いのではないかと思います。

[0.001, 0.01, 0.1, 1, 10, 100]

 ここで、0.001が出る頻度に対して0.002はそこそこ出てきてほしいけど、100の頻度に対して100.001は出てきてほしくなくて200とかになってほしい訳です。要するに、上の方になるほど出てくる確率が下がる分布がほしい訳で、こういうのは指数分布が向いている気がします。ということを踏まえて、scipy.stats.exponを使うことにします。

scipy.stats.expon — SciPy v1.2.1 Reference Guide

 scipyのこの辺のドキュメントの説明はお世辞にもわかりやすいとは言えないんですが、指数分布のパラメータは基本的に \lambda一つだけであり、ドキュメントの記述によると

A common parameterization for expon is in terms of the rate parameter lambda, such that pdf = lambda * exp(-lambda * x). This parameterization corresponds to using scale = 1 / lambda.

 ということらしいので、これに従って \lambdaを設定します。指数分布の期待値は \frac{1}{\lambda}なので、要するに平均にしたいあたりをscaleに指定すればいいことになります。

 ちょっと確認してみます。

>>> import matplotlib.pyplot as plt
>>> result = stats.expon.rvs(scale=0.1, size=1000)
>>> plt.figure()
<matplotlib.figure.Figure object at 0x7fe9e5e4b6a0>
>>> plt.hist(result)
(array([639., 222.,  76.,  44.,  11.,   3.,   3.,   1.,   0.,   1.]), array([7.33160559e-07, 9.56411007e-02, 1.91281468e-01, 2.86921836e-01,
       3.82562203e-01, 4.78202571e-01, 5.73842938e-01, 6.69483306e-01,
       7.65123673e-01, 8.60764041e-01, 9.56404408e-01]), <a list of 10 Patch objects>)
>>> plt.savefig("result1.png")
>>> result = stats.expon.rvs(scale=1, size=1000)
>>> plt.figure()
<matplotlib.figure.Figure object at 0x7fe9e5e1bba8>
>>> plt.hist(result)
(array([490., 246., 133.,  66.,  30.,  23.,   2.,   6.,   3.,   1.]), array([2.21631197e-04, 6.36283260e-01, 1.27234489e+00, 1.90840652e+00,
       2.54446814e+00, 3.18052977e+00, 3.81659140e+00, 4.45265303e+00,
       5.08871466e+00, 5.72477629e+00, 6.36083791e+00]), <a list of 10 Patch objects>)
>>> plt.savefig("result2.png")

result1.png
result1.png
result2.png
result2.png

 上手く行っているようですが、指数分布の場合は分散が \frac{1}{\lambda ^2}で要するに期待値の2乗なので、場合によってはいわゆる中央値より少し大きめにしたり、逆に小さめにしたいと思うこともあるかもしれません。

 というあたりを踏まえて、いよいよRandomizedSearchCVでやってみます。基本的な使い方はコードを見れば分かる通りで、paramsの値にscipyの確率分布を渡すこと、n_iterで試す回数を指定できること以外大きな違いはありません。n_iterですが、今回は30回やってみます。

コード

import time

import pandas as pd

from scipy import stats
from sklearn.svm import SVC
from sklearn.model_selection import RandomizedSearchCV
from sklearn.datasets import load_digits

def main():
    digits = load_digits()
    svm = SVC()
    params = {"C":stats.expon(scale=1), 
              "gamma":stats.expon(scale=0.01)}
    clf = RandomizedSearchCV(svm, params, cv=5, iid=True,
                             return_train_score=False, n_iter=30)
    t1 = time.time()
    clf.fit(digits.data, digits.target)
    t2 = time.time()
    print("{:.2f}秒かかった".format(t2 - t1))
    result_df = pd.DataFrame(clf.cv_results_)
    result_df.sort_values(
        by="rank_test_score", inplace=True)
    print(result_df[["rank_test_score", 
                     "param_C",
                     "param_gamma",
                     "mean_test_score"]])

if __name__ == "__main__":
    main()

結果

63.04秒かかった
    rank_test_score     param_C  param_gamma  mean_test_score
23                1     2.62231  0.000465122         0.972732
15                2     1.85208   0.00195101         0.967168
14                3    0.716893   0.00232262         0.963829
26                4    0.633723  0.000583885         0.960490
20                5    0.534343   0.00280065         0.952699
24                6     1.33912   0.00344179         0.949360
2                 7     4.06141   0.00355782         0.947691
21                8    0.460892   0.00307683         0.945465
19                9     1.78068   0.00481671         0.912632
10               10     3.26205   0.00679844         0.817474
29               11    0.288305   0.00525321         0.738453
1                12     1.01437   0.00874084         0.731219
3                13    0.385384   0.00590435         0.725097
5                14     0.61161   0.00702304         0.701169
12               15    0.317451   0.00597252         0.626600
6                16    0.122253   0.00444626         0.521425
4                17      1.1791    0.0187758         0.373400
28               18    0.960299    0.0174957         0.277685
7                19     1.91593    0.0219393         0.272120
18               20     1.17152    0.0264625         0.209794
13               21     1.48332    0.0297287         0.180301
22               22    0.496221    0.0110038         0.178631
11               23    0.747978    0.0178759         0.130217
17               24    0.549571    0.0147944         0.123539
0                25   0.0357021   0.00627669         0.116305
25               26   0.0551073   0.00739485         0.114079
27               27   0.0852409   0.00982068         0.111297
9                27    0.148459    0.0101576         0.111297
16               29    0.818395    0.0521121         0.102393
8                29  0.00447628    0.0442941         0.102393

 最良スコアそのものは実はGridSearchCVと(たまたま)同じですし、時間がケチれるということもないのですが、見ての通りある程度アタリのついている状態でうまい分布を指定してあげれば「最良パラメータの周囲を細かく探索する」といったカスタマイズが可能になります。

 とりあえずこれでどの辺りがいいのかは大体わかったので、今度はその辺りを平均にしてやってみます。2行、次のように書き換えただけです。

改変箇所

    params = {"C":stats.expon(scale=3), 
              "gamma":stats.expon(scale=0.001)}

結果

22.06秒かかった
    rank_test_score      param_C  param_gamma  mean_test_score
28                1      4.17858  0.000512401         0.974402
3                 2      3.64468  0.000429315         0.973845
14                3      6.71633  0.000842716         0.973289
22                3      2.69443   0.00116543         0.973289
6                 3      5.40707  0.000698026         0.973289
13                3      2.80969   0.00118559         0.973289
1                 7      3.00941   0.00135159         0.972732
11                7      2.10982  0.000993251         0.972732
12                7      3.73687   0.00112881         0.972732
2                10      1.18243   0.00144449         0.971063
24               11      1.86039   0.00167388         0.970506
16               11      1.25451   0.00150173         0.970506
10               13      2.45206  0.000338623         0.968837
0                14     0.827324   0.00171363         0.967724
25               15     0.651184   0.00175566         0.965498
29               16      1.24774  0.000424962         0.964385
20               17     0.625333  0.000705793         0.963829
15               18      7.14503  0.000158783         0.962716
27               19     0.461709   0.00106818         0.962159
5                20      2.82219   0.00280609         0.961046
4                21      1.47712  0.000229419         0.959377
7                22     0.536423    0.0005259         0.958264
8                23      2.39408  9.90463e-05         0.953812
18               24      3.53976  6.63105e-05         0.953255
21               24      3.46599  6.85291e-05         0.953255
23               26      3.43756  5.60631e-05         0.952142
26               27       12.188  2.79396e-05         0.951586
9                28     0.331784   0.00038362         0.951029
17               29      0.73821   4.9347e-05         0.930440
19               30  0.000110158  0.000966082         0.141347

 今度はGridSearchCVのときより良いスコアになるパラメータがいくつか出ました。まあ、digitsの大して多くもないデータ数で0.001の差を云々しても「誤差じゃね?」って感じですが・・・

 二段グリッドサーチでやっても同じことはできますが、グリッドを手打ちで入力したりといった手間がかかります。その点、RandomizedSearchCVは分布のパラメータを打ち込めば勝手にやってくれるので、大変助かるものがあります。

結論

 使えます。特に、グリッドサーチのグリッドを手打ちで入れるのが面倒くさいという人に向いています。ただし、どういう分布が良いのかは知識として持っていないといけません。まあ、わからなければ一様分布でアタリをつけて正規分布にして・・・とかでもなんとかなるでしょう。

 ただし、それなりに工夫する要素があるので、玄人向きです。kaggleで使われたりするのも納得がいきます。また、確率である以上ある程度の回数は回さないと安定しないので、そのへんにも注意した方が良いと思います(それでも数が増えてくるとグリッドサーチできめ細かくやるより全然マシですが・・・)。

numpyやpandasでThe truth value of ... is ambiguous.のようなエラーが出たときの対処

概要

 条件式を使って生成したようなboolのnumpy配列を使っていると、次のようなエラーが出ることがあります。

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

 また、pandasのSeriesやDataFrameでも同様のエラーが発生する場合があります。

ValueError: The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all().
ValueError: The truth value of a DataFrame is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all().

 この記事では、こういったエラーの原因と対処法について説明します。

原因

 numpyやpandasの複数の値を持つコレクション型は、そのままbool値に変換できないことになっています。なぜか? というと、そういったケースではコレクション全体の真理値は曖昧(ambiguous)だからです。

 ValueErrorのメッセージでは以下のような変換が提案されています。

  • a.any()かa.all()を使う

 anyは要素すべてのOR、allは要素すべてのANDです。このような方法で変換するのが適切であれば、そうしてくれということです。

  • a.empty(pandasのオブジェクトのみ)

 空であるかどうか。

  • a.bool()(pandasのオブジェクトのみ)

 ブール型の単一の要素を持つ場合、その真理値が返ります。それ以外のケースではValueErrorになるようです。

  • a.item()(pandasの)

 要素数1のSeriesに対してその唯一の要素を返します。

 なんとなく雰囲気がつかめてきましたね。何らかの方法で1つにまとめるなり、1つだけ取り出すなりしてくれと言っている訳です。でもまあ、この通りにすればよいかというと必ずしもそうではなく、ケースバイケースで対処する必要があります。

 余談ですが、このような挙動をするのは私が確認した範囲ではnumpyやpandasのコレクション型のみで、組み込みのlist型などではこの挙動にはなりません。どうなるかというと、皆さんよくご存知の標準的なpythonと同等の真理値判定が行われます。

4. 組み込み型 — Python 3.6.8 ドキュメント

 numpyやpandasの実装は、大雑把に言えば__bool__が例外を送出するようにしてあるというだけのものです。たとえばpandasはこんな感じ↓です。

https://github.com/pandas-dev/pandas/blob/master/pandas/core/generic.py#L1498

 どうしてわざわざこんなことをしているのか? というと、素のpythonと同等の判定だと(たとえば空だとFalseとか)中身のデータに対する真理値だと思いこんで処理する人が出てきたときにバグが発生するので、配慮して__bool__を潰してあるのだと思われます。要するに、こんな親切なメッセージまで出してくれるのは「優しさ」なので、「なんだこの腹立つエラーは」とか思ってはいけません。

対処

 対処法ははっきり言ってケースバイケースです。このエラーはいろいろな原因で発生するので、状況にあった対処をする必要があります。幾つか紹介してみます。

if文の条件式にそのまま書いた

 if文は内部で条件式のbool値への変換を行います。ということは、このエラーが発生する可能性があります。

if numpyやpandasのコレクション型:
    ....

 これに関しては、はっきり言ってこういうことをやろうとするのが悪いです。何かしら考えが間違っていると思います。要素を一つ一つループで取り出すべきかもしれませんし、indexingで処理するべき状況の可能性もあります。あるいはそもそも期待と違うものが変数に代入されている可能性もあります。全般的にデバッグしてください。

andやorを使った

 andやorなどの演算子は「pythonのbool」に対して働くので、numpy配列、pandasのSeriesやDataFrameに対して使うのは基本的に不適当です。これらは内部の判定時にboolへの変換を行うので、この記事で取り上げているエラーを発生させる要因になります。

print(np.array([True, False]) and np.array([True, True]))
# => ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

 この場合、期待しているものは&演算子と|演算子で得られるかと思われます。

print(np.array([True, False]) & np.array([True, True]))
# => [ True False]

複数の条件式の組み合わせで発生した

 上の2つは割と初歩的ですが、このケースは割と上級者の人でもハマるときはハマります。

 上述の通り、配列同士の論理演算は&と|でできるのですが、これは==のような比較演算子より優先順位が強い演算子になります。まあ、本来はbit演算用の演算子なのを流用しているので仕方ないのですが、結果的に評価の順序が狂ってエラーになるケースがあります。

6. 式 (expression) — Python 3.7.3rc1 ドキュメント

a = np.array([0,1,2])
b = np.array([3,4,5])
print(a == 1 | b == 3)
# => ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

 これはこう書いたのと同じです。

print(a == (1 | b) and (1 | b) == 3)

 pythonの比較演算子の仕様がちょっと特殊で、連結するとandで連結したのと同じになることに留意。

6. 式 (expression) — Python 3.7.3rc1 ドキュメント

 幸いかっこを付けて比較的容易に優先順位を制御できます。

a = np.array([0,1,2])
b = np.array([3,4,5])
print((a == 1) | (b == 3)) # => [ True  True False]

 しかしまあ、ちょっと困った話です。けっきょくpythonの世界であれこれ処理しようとするとこういう細かい齟齬が出てくるので、pandasにqueryメソッドがあったりするのもむべなるかなという感があります。かといってqueryの書き方を覚えるのも面倒くさいんですが。

まとめ

 割とよく見かけるエラーですが、基本的に何かがおかしくなってしまったときに出るもので、ケースバイケースで対応することになるからかあまり掘り下げた解説がなかったので、書いてみました。この記事で深く掘り下げられているかというと少し微妙な感がありますが・・・

 とにかく「複数の値が入っているから真理値が一意に決められない(ことになっている)」というのがミソなので、それだけは覚えておきましょう。

ファイルオブジェクトのcloseはflushも行う。確実にしたければfsync

 以前、「ファイルオブジェクトのcloseメソッドは同時にflushも行う」ことを知りました。

 どうやらcloseするときは内部でflushメソッドが呼ばれるようです。

このストリームをフラッシュして閉じます。
io --- ストリームを扱うコアツール — Python 3.7.3rc1 ドキュメント

 ということは、「ファイルを閉じた後即座に反映させたいので、flushメソッドを呼んでおく」という操作は必要ないようです。closeすれば同時にflushされるので。

 余談ですが、ファイルオブジェクトのflushを呼べばそれで確実に処理が終わったタイミングでディスクに書かれているか? というと、OS側でバッファリングされるので、保証はされないということらしいです。

operating system - does close() imply flush() in Python? - Stack Overflow

 flushはあくまでも「pythonプロセスからOSにちゃんと渡した」というところまでしか保証しないのでしょうか・・・。

 対策はfsyncと呼ばれるシステムコールを使うことです。pythonでもosモジュールにインターフェースがあります。

参考:ファイル出力したそのデータ、本当にディスクに書き出されてますか? - らっちゃいブログ

Python の ファイルオブジェクト f を使う場合、f の内部バッファを確実にディスクに書き込むために、まず f.flush() を、その後 os.fsync(f.fileno()) を実行してください。
os --- 雑多なオペレーティングシステムインタフェース — Python 3.7.3rc1 ドキュメント

 あ、これはflushしないと駄目なんですね。先にcloseしないので当然か。

import os

with open("sample.txt", "w") as f:
    f.write("hoge")
    f.flush()
    os.fsync(f.fileno())

 このようにすると処理を抜けた時点でディスクに書かれていることが保証されると思います。また、メタデータ更新を行わない(ドキュメントには「強制しない」と書いてあるが、実際にどういう挙動なのかは謎)fdatasyncという類似関数も用意されています。上に載せた参考リンクで紹介されています。

scipy.optimize.curve_fitを使っていろいろな関数にフィットさせてみる

はじめに

 scipy.optimize.curve_fitを使うと曲線あてはめができます。いろいろな関数にフィッティングさせてみて、うまくいくかどうか試してみます。

scipy.optimize.curve_fit — SciPy v1.2.1 Reference Guide

f(x) = x + a

 ただの足し算。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return x + a

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, a in enumerate([2, 5, 14]):
        y = func(x, a) + np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label="a={:.4f}".format(*params))
    plt.legend()
    plt.savefig("result1.png")

if __name__ == "__main__":
    main()

result1.png
result1.png

 問題なさそうです。

f(x) = a * x

 掛け算。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return x * a

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, a in enumerate([1, 1.5, 2]):
        y = func(x, a) + np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label="a={:.4f}".format(*params))
    plt.legend()
    plt.savefig("result2.png")

if __name__ == "__main__":
    main()

result2.png
result2.png

 これも問題なし。

f(x) = a*x + b

 一次式です。これくらいはできないと困ります。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a, b):
    return x * a + b

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, pt in enumerate([[1, 2], [3, 4], [5,6]]):
        y = func(x, *pt) + np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label="a={:.4f}, b={:.4f}".format(*params))
    plt.legend()
    plt.savefig("result3.png")

if __name__ == "__main__":
    main()

result3.png
result3.png

 難なくこなしました。

4次多項式

 次数の少し大きめの多項式でやってみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a, b, c, d, e,):
    return a*x**4 + b*x**3 + c*x**2 + d*x + e

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, pt in enumerate([[1,2,3,4,5],
                            [2,3,4,5,6],
                            [3,4,5,6,7]]):
        # 値のスケールが大きいのでノイズの大きさを調整
        y = func(x, *pt) + 1000*np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label=("a={:.3f} b={:.3f}"
                        "c={:.3f} d={:.3f}"
                        "e={:.3f}").format(*params))
    plt.legend()
    plt.savefig("result4.png")

if __name__ == "__main__":
    main()

result4.png
result4.png

 当たり前ですが、全体の傾向は高次の係数に支配されるため、低次の係数ほど元の値と乖離します。それでもノイズが小さければ使えるんですが。

aのx乗

 なんとなく難しそうな気がしました。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return a**x

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, pt in enumerate([[1], [1.2], [1.4]]):
        y = func(x, *pt) + np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label=("a={:.3f}").format(*params))
    plt.legend()
    plt.savefig("result5.png")

if __name__ == "__main__":
    main()

result5.png
result5.png

 意外となんとかなりました。けっこういい精度のように見えます。

xのa乗

 こちらもやってみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return x**a

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, pt in enumerate([[1], [1.2], [1.4]]):
        y = func(x, *pt) + np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label=("a={:.3f}").format(*params))
    plt.legend()
    plt.savefig("result6.png")

if __name__ == "__main__":
    main()

 上のコードを実行しようとすると、エラーを吐きます。

***.py:6: RuntimeWarning: invalid value encountered in power
  return x**a
/***/site-packages/scipy/optimize/minpack.py:715: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)
Traceback (most recent call last):
  File "***.py", line 22, in <module>
    main()
  File "***.py", line 14, in main
    params, cov = curve_fit(func, x, y)
  File "/***/site-packages/scipy/optimize/minpack.py", line 654, in curve_fit
    ydata = np.asarray_chkfinite(ydata)
  File "/***/site-packages/numpy/lib/function_base.py", line 1033, in asarray_chkfinite
    "array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs

 区間が駄目なのだろうと思って、0を含まない区間に変えます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return x**a

def main():
    n = 50
    x = np.linspace(1, 10, n)
    plt.figure()
    for i, pt in enumerate([[1], [1.2], [1.4]]):
        y = func(x, *pt) + np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label=("a={:.3f}").format(*params))
    plt.legend()
    plt.savefig("result6.png")

if __name__ == "__main__":
    main()

result6.png
result6.png

 今度はできました。0を何乗しても0なので別に問題ない気がしますが、内部の計算に失敗するのかもしれません。もしかして微分して0になる系は厳しい?

対数の底

 対数の底をパラメータで探ってみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return np.log(x)/np.log(a)

def main():
    n = 50
    x = np.linspace(1, 10, n)
    plt.figure()
    for i, pt in enumerate([[2], [3], [4]]):
        y = func(x, *pt) + 0.1*np.random.randn(n)
        params, cov = curve_fit(func, x, y)
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label=("a={:.3f}").format(*params))
    plt.legend()
    plt.savefig("result7.png")

if __name__ == "__main__":
    main()
***.py:6: RuntimeWarning: divide by zero encountered in true_divide
  return np.log(x)/np.log(a)
***.py:6: RuntimeWarning: invalid value encountered in true_divide
  return np.log(x)/np.log(a)
/***/site-packages/scipy/optimize/minpack.py:787: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

 たぶん初期値が駄目なパターンなので、初期値として1.1を与えてみます。それ以外は変更なし。

        params, cov = curve_fit(func, x, y, p0=[1.1])

result7.png
result7.png

 問題なくフィッティングできました。

 p0を指定しなかった場合、パラメータの初期値は0になります。関数によっては0から始まると困ったことになるので、適当な初期値を与えることが重要です。

シグモイド関数

 割と複雑な関数。シグモイドで曲線あてはめしたいというのはそれなりに実用的なユースケースでしょう。

 yにノイズを足すと0より小さかったり1より大きい値が出てきてしまうので、xをずらしました。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a):
    return (np.tanh(a*x/2)+1)/2

def main():
    n = 50
    x = np.linspace(-10, 10, n)
    plt.figure()
    for i, pt in enumerate([[1], [10], [100]]):
        y = func(x + 0.2*np.random.randn(n), *pt)
        params, cov = curve_fit(func, x, y, p0=[1.1])
        plt.scatter(x, y, c="rgb"[i])
        plt.plot(x, func(x, *params), c="rgb"[i],
                 label=("a={:.3f}").format(*params))
    plt.legend()
    plt.savefig("result8.png")

if __name__ == "__main__":
    main()

 このまま実行すると、数回に一回くらいうまくフィットしないことがあります。

/*/lib/python3.5/site-packages/scipy/optimize/minpack.py:787: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

 そのときの結果がこんな感じです。

result8.png
result8.png

 データが悪い。初期値を真値に近づけると多少改善するかもしれません。

まだやっていないもの

  • 各種確率密度関数など(正規分布、指数分布、多項分布あたり)
  • シグモイド関数に似ているもの(ゴンペルツなど)

 このあたりは現時点ではやっていませんが、基本的には同様にできると考えられるので、気が向いたら追記します。

【python】numpyで任意の底でlog

 任意の底でlogを計算したいときがあります。

 結論から言うと、そういう関数は用意されていません。log, log10, log2はあるんですが。

Mathematical functions — NumPy v1.16 Manual

 ほしければ、自分でlog(底)で割ってあげます。

>>> import numpy as np
>>> a = 59**5
>>> np.log(a)/np.log(59)
5.0

 対数の性質を思い出せばそれほど難しくないのですが、デフォルトであると便利ですね。math.logなら一発でできるんですが。

>>> import math
>>> math.log(a, 59)
5.0

 おまけ。

 最初、すごく大きな値でやろうとしたらエラーが出ました。

>>> a = 59**100
>>> np.log(a)/np.log(59)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'int' object has no attribute 'log'
>>> a
1216749465718417874559491998513561121966072058691773438388982975785993854256756157494173512672915757033875049659241870497099358702416318491812748479721683887882543089775014614001

 これは多倍長整数だから表現できているのであって、int64とかでは表現できないためエラーになるということですね。floatに直したら行けました。

>>> float(a)
1.2167494657184178e+177
>>> np.log(float(a))/np.log(59)
100.0

 精度は有効桁数に縛られますが(doubleの有効桁数はたかだか15桁)、割とどうでも良い。

【python】scipyで線形最小二乗法

概要

 scipyのscipy.optimize.lsq_linearで最小二乗法が使えます。

使い方

 最低限必要な引数は、

  • A

 いわゆる説明変数です。基本的には(データ数, 次元数)のshapeでいいのですが、バイアス項を入れたければすべて1にした列が要ります。

  • b

 いわゆる目的変数です。(データ数,)のshapeで入れます。

 の2つだけです。他に色々なパラメータを設定できます。

 返り値はscipy.optimize.optimize.OptimizeResultという型で返り、表示すると最適化の結果をレポートしてくれます。単に係数を知りたければxというフィールドを見れば良いです。

scipy.optimize.lsq_linear — SciPy v1.2.1 Reference Guide

使ってみる

 簡単に試してみます。必要に応じてコメントを入れます。

>>> import numpy as np
>>> X = np.random.rand(100, 3)  # ランダムなXを作成
>>> X[:,2] = 1  # バイアス項を設定
# 適当なyを作成。正規分布するノイズをまぶす
>>> y = 2*X[:,0] + 3*X[:,1] + 4 +0.5*np.random.randn(100)  
>>> from scipy import optimize
>>> result = optimize.lsq_linear(X, y)
>>> result
 active_mask: array([ 0.,  0.,  0.])
        cost: 12.474073406201914  # コスト関数の値
        # 残差
         fun: array([-0.62398486,  0.3072979 , -0.10238148, -0.4993848 ,  0.0734566 ,
       -0.06842613, -0.78150259, -0.28869525, -0.29234346, -0.09920162,
        0.67799186,  0.47657423, -0.14040916,  0.22752255,  0.35423127,
        0.90124369,  0.31145881,  0.12874122,  1.03422167, -0.30381206,
       -0.22739989, -0.02065523,  0.3710351 ,  0.12505248,  0.29515137,
       -0.20530665, -0.19887154,  0.2782884 , -0.03371473, -0.82749514,
        0.4542083 ,  0.00581736, -0.11176233,  0.7640272 ,  0.59850616,
        0.21529376,  0.29703533,  0.85004091,  0.73353464, -0.07128385,
       -0.86266412, -0.02184265, -0.57813797,  0.58658777, -0.11631637,
        0.38530319, -1.52933752, -0.2369261 , -0.30857852,  0.68003027,
        0.07390349, -0.54587462,  0.22403256, -0.08884161,  0.40610932,
       -0.41575813,  0.14871697,  0.34747144, -0.54337788, -0.02348569,
        0.44514965, -0.03702159,  0.07137893, -0.55835038,  0.65997585,
        0.42452063,  0.71298796, -0.019805  ,  0.32504799,  0.51807313,
        0.09903619, -0.99172989, -0.04408202, -0.20914935,  0.37428201,
       -0.49526746, -0.87998355, -1.03734365, -0.44080683, -0.5162824 ,
        1.19186206, -0.64073783, -0.10657965, -0.07972209, -0.36248635,
        0.47518401, -0.12845207, -0.44213326, -0.67951796,  0.44330662,
       -0.25102293,  0.5503437 ,  0.10920638,  0.21309121, -0.06616317,
       -0.10048913, -0.32081914, -0.1639193 , -1.00984118,  0.803146  ])
     message: 'The unconstrained solution is optimal.'
         nit: 0
  optimality: 5.5777604757167865e-13
      status: 3
     success: True 
          # 係数
           x: array([ 2.03207968,  3.19444993,  3.90466908])
>>> result.x
array([ 2.03207968,  3.19444993,  3.90466908])
|<