2020-08-24

pythonで標準出力を取得する。

f = io.StringIO()
with redirect_stdout(f):
    help(pow)
s = f.getvalue()

公式ドキュメント

終わり！

とするのは味気ないのでなぜこんなことを調べたのか補足。

生存時間解析のライブラリの[lifelines] (https://lifelines.readthedocs.io/en/latest/index.html)を使っているたのですが、fittingした結果をprint()で返す仕様になっていました。
jupyter notebookで使う間はこの仕様でも問題なく動作するのですが、streamlitへ出力する場合には、標準出力は具合が悪いです。
例えばlifelineのCox比例ハザードモデルを実行したとします。

from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi

rossi = load_rossi()

cph = CoxPHFitter()
cph.fit(rossi, duration_col='week', event_col='arrest')

cph.print_summary()

このときの print_summary() 関数はソースコードを見るとprint() 関数で標準出力しています。

なので、 redirect_stdout() を使って、

from io import StringIO
from contextlib import redirect_stdout


f = StringIO()
with redirect_stdout(f):
    cph.print_summary(style='html')

とすれば、htmlで得られるので、あとはstreamlitを使って、

import streamlit as st

st.markdown(f.getvalue(), unsafe_allow_html=True)

として出力できます。streamlitのissueを見ていると、

github.com

というissueもあるようなので、そのうちstreamlit側で解決されるかも。

2020-08-02

pythonの相対インポートがわからなくなったメモ

たまにわからなくなるので、メモ書きです。

pytorchのgithubを読んでいると、次のような相対importが出てきます。

from .module import Module

https://github.com/pytorch/pytorch/blob/master/torch/nn/modules/conv.pyのLine 10とか。

これを参考に、自分でも相対インポートをしようと思いました。

src/  
├── bar.py  
└── foo.py

foo.py

def foo():
    print('foo')

bar.py

from .foo import foo

def bar():
   print('bar')
   foo()

if __name__ == '__main__':
     bar()

python bar.py を実行しようとしたのですが、このときには

Traceback (most recent call last):  
  File "bar.py", line 1, in <module>  
    from .foo import foo  
ImportError: attempted relative import with no known parent package

とエラーが出ます。
非常に初歩的な内容なのですが、.を用いた相対インポートを行う場合には、親ディレクトリが必要なようです。
なので、

src
├── main.py
└── module`
    ├── __init__.py
    ├── bar.py
    └── foo.py

のようにディレクトリ構成を変更して、

module1/__init__.py

from .bar import bar
from .foo import foo

__all__ = [foo, bar]

main.py

from module1 import foo, bar

if __name__ == '__main__':
    foo()
    bar()

として、main.pyを実行すれば、

foo
bar
foo

と返り値が得られました。

なので、相対インポートを使いたいときには親ディレクトリを作って使いましょう。

実行ファイルと同じ並びに相対インポートしたいpythonファイルがある場合には、. を使わず、 bar.py

from foo import foo

def bar():
   print('bar')
   foo()

if __name__ == '__main__':
     bar()

とすれば、インポートできます。

2019-09-30

pytorchでpretrained modelの最終層だけlrを変更する。

torchvisionで呼び出した最終層だけlrを変えるやり方がわからなかったので調べてみました。 pytorchの公式ドキュメントに記載があるので，そちらを参照して頂いたほうが正しい情報になると思います。

まず適当なmodelを呼び出します。

from torchvision import models

model = models.resnet18(pretrained=False)

lrを変えるのは，optimizerの部分で定義します。公式ドキュメントをフォローしています。

import torch.optim as optim
optimizer = optim.SGD([
                                        {'params': model.parameters()},
                                        {'params': model.fc.parameters(), 'lr': 1e-3}
                                       ], lr=1e-2, momentum=0.9)

で行けるようです。あとは学習させるだけです。

2019-09-28

PILで取り込んだ画像をnumpyに変換してmaskと重ねる

普段，私は画像はPIL(Pillow)で取り込んでいます。色々あって，semantic segmentationを取り組んでいるのですが，やっぱりmaskデータと重ねて描写してみたくなります。いまいちやり方がわからなかったので，自分なりにググってやってみたことを書いておきます。単なる備忘録です。絶対もっといいやり方があると思います。

参考にしたのは，

stackoverflow.com

です。

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt

#データ取り込み
image = Image.open(img_data)
mask = np.array(mask_data)

# numpyへの変換
image = np.asarray(image)

# figureサイズ操作
plt.figure(figsize=(XX, XX), dpi=XX)

# 描写
plt.imshow(image)
plt.imshow(mask, "jet", alpha=0.2)

です。最後の引数"jet"は色の好みで，alphaで透明度を変更できます。

2019-05-01

ベイズ線形回帰

python

動機

MCMCで計算をすると，カルマンフィルタを実装したときは実感できたベイズ更新の部分がよくわからないという僕の気持ちの問題がありました。
そこで，須山さんの「ベイズ推論による機械学習入門」

機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)

作者: 須山敦志,杉山将
出版社/メーカー: 講談社
発売日: 2017/10/21
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

の3章で，学習と推論を解析的に導出していたため，自分で計算・実装を行い理解を深めたいというのがこの記事を書く動機です。

線形回帰の式

１次関数での線形回帰

想像しやすいように，１次関数での線形回帰を行います。式は $y = a x + b+ \varepsilon_r$ です。
$\varepsilon_r$ はノイズで，回帰分析の事後分布として求めたい値は， $a, b$ となります。
ノイズは簡単のために正規分布に従うと仮定して，
$\varepsilon_r \sim N(\varepsilon_r | 0, \lambda^{-1})$
と定義します。テキストに習って $\lambda$ は既知であるとします。
ある一つの観測値 $y_k$ とある一つの入力値 $x_k$ が得られたときに，事後分布 $p(a| y_k, x_k), p(b|y_k, x_k)$ が導出できれば，学習と予測ができることになります。

まず，簡単のために，bは既知であるとして，aの事後分布を解析的に求めてみます。
ベイズの定理より，求めたい事後分布 $p(a| y_k, x_k)$ は，
$p(a| y_k, x_k) = \frac{p(y_k|x_k, a) p(a)}{p(y_k, x_k)}$
と書けます。
このときに，事前分布 $p(a)$ は正規分布していると仮定し
$p(a) = N(a | m_a, \Lambda_a^{-1})= \frac{1}{\sqrt{2\pi\Lambda_a^{-1}}}\exp\big( -\frac{1}{2}\Lambda_a(a-m_a)^{2} \big)$
と書けるとします。ここで， $m_a, \Lambda_a$ は事前分布の平均と分散で固定値で与えるハイパーパラメータです。

さて，事前分布 $p(a)$ と入力値 $x_k$ が与えられた元で，観測値 $y_k$ が得られる度合いを表す尤度を考えます。
先程，ノイズが正規分布に従うとしたため，
$p(y_k|x_k, a) = N(y_k | a * x_k +b, \lambda^{-1} )$
と書くことができます。ここで，この正規分布を真面目に書いてみると
$N(y_k | a * x_k+b, \lambda^{-1}) = \frac{1}{\sqrt{2\pi\lambda^{-1}}}\exp\big( -\frac{1}{2}\lambda(y_k-a*x_k-b)^{2} \big)$
となります。
今は $a$ の事後分布だけが気になっているので，計算を簡単にするために事前分布，尤度をそれぞれの対数( $\log$ )を取って， $a$ について整理していきます。

事前分布

$\log p(a) = \log \frac{1}{\sqrt{2\pi\Lambda_a^{-1}}}\exp\big( -\frac{1}{2}\Lambda_a(a-m_a)^{2} \big)$
$= -\frac{1}{2} \log (2\pi\Lambda_a^{-1} ) -\frac{1}{2}\Lambda_a(a-m_a)^{2}$

第一項は $a$ に関係のない項なので，第二項を展開して， $a$ に関する項を昇順に並べていきます。
$-\frac{1}{2}\Lambda_a(a-m_a)^{2} = -\frac{1}{2} \Lambda_a a^{2} + \Lambda_a m_a a -\frac{1}{2} \Lambda_a m_a^{2}$ (1)
となります。

尤度

次に，先程定義した尤度も同様に計算することで，
$-\frac{1}{2} \lambda_a x_k ^{2} a^{2} + \lambda (y_k -b) x_k a -\frac{1}{2} \lambda (y_k -b )^{2}$ (2)
と書けます。

事前分布と尤度の積

ベイズの定理の分子である $p(y_k|x_k, a) p(a)$ は事前分布と尤度の積ですので，対数を取ったあとは和になります。
なので，(1)と(2)の和を計算して， $a$ について昇順に並べることで，事後分布の $a$ についてはどのような形になるのか想像することができます。和を取ると $-\frac{1}{2} ( \Lambda_a + \lambda_a x_k ^{2} ) a^{2} +( \Lambda_a m_a + \lambda (y_k -b) x_k )a -\frac{1}{2} ( \Lambda_a m_a^{2} + \lambda (y_k -b )^{2})$ (3)
と書けます。
この形は事前分布や尤度で仮定した正規分布と全く同じ形になっていることから，事後分布も正規分布であると考えることができます。*1
分母は，周辺尤度を求めれば，正規化されていると思います。
このあたりの数式は天下り的に計算はしませんが，任意の事後分布の正規化の数式などは例えば，

ノンパラメトリックベイズ点過程と統計的機械学習の数理 (機械学習プロフェッショナルシリーズ)

作者: 佐藤一誠
出版社/メーカー: 講談社
発売日: 2016/04/20
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

こちらの本などに書かれているように思えます。
さて，事後分布も正規分布であるとして，一点を観測したあとの傾き $a$ の事後分布は
$p(a| y_k, x_k) = N(a| ( \Lambda_a + \lambda_a x_k ^{2} ), \frac{\Lambda_a m_a + \lambda (y_k -b) x_k}{( \Lambda_a + \lambda_a x_k ^{2} )})$ (4)
として計算できます。*2

今，観測点が一点だけの場合を考えていましたが，複数観測する場合には，式(3)は更に一般化することができて， $-\frac{1}{2} ( \Lambda_a +\sum_k \lambda_a x_k ^{2} ) a^{2} +( \Lambda_a m_a + \sum_k \lambda (y_k -b) x_k )a -\frac{1}{2} ( \Lambda_a m_a^{2} + \sum_k \lambda (y_k -b )^{2})$ (3)
と，書けます。このとき，(4)は
$p(a| y_k, x_k) = N(a| \frac{\Lambda_a m_a + \sum_k \lambda (y_k -b) x_k}{( \Lambda_a + \sum_k \lambda_a x_k ^{2} )}, ( \Lambda_a + \sum_k \lambda_a x_k ^{2} ))$ (5)
となりますので，この式を使うことで，一点一点を更新していく様子がみることができます。

実装

実装コードはこちらにあります。

github.com

傾き $a=5$ , 切片 $b=3$ , ノイズが $\sigma=1$ の正規分布に従って生成されるとして，データをランダムに200点ほど作成します。

f:id:Leo0523:20190501223555p:plain — 図1: 観測データ

図1のような，データが得られます。

データが得られたので，上の計算で行ったように，切片aをベイズ線形回帰で計算してみます。
事前分布は，中心ピークを4.5，標準偏差 $\sigma=2$ の正規分布を仮定します。

f:id:Leo0523:20190501223410p:plain — 図2: 事前分布

この仮定の元，上記の計算を実行して，事後分布が毎回どのような形状になるのか計算してみます。
結果は次の動画で，

f:id:Leo0523:20190501224331g:plain — ベイズ更新

一点ごとの更新となりますので，徐々に予測の幅が狭くなり，値が真値である5に収束していく様子が見えます。
これは僕が当初見たかったベイズ線形回帰におけるベイズ更新の様子になります。

現在の仮定では事前分布の幅が広く，ピークも真値に近いので，パラメータを変えて実験を行ってみました。

事前分布：ピーク3, 標準偏差 $\sigma=0.5$ , 観測データ点数100点

f:id:Leo0523:20190501225109g:plain — ベイズ更新2

事前分布：ピーク3, 標準偏差 $\sigma=0.1$ , 観測データ点数100点

f:id:Leo0523:20190501225247g:plain — ベイズ更新3

標準偏差が小さいときには，事前分布に引っ張られて真の値になかなかたどり着かない様子がこれらの動画から見ることができました。
事前分布，大事！

まとめ

教科書に従って事後分布を解析的に導出した。
導出した式を使ってベイズ更新の様子を描写した。
（数式的には自明だが）事前分布の標準偏差が小さいほど真値を予測するには観測データが沢山必要であることが動画からわかった。

*1:ほんまかいな？数学的な証明はわかりません

*2:詳細な計算は須山さんのベイズ推論による機械学習入門を読んでください

2019-01-13

StanとRでベイズ統計モデリング,　４章の計算をしてみた

PyStan

概要

StanとRでベイズ統計モデリングの4章のベイズ信頼区間と予測区間をpythonを自分で実装してみました。

github.com

実験したコードはこちらです。

動機

ベイズ統計を実務で使うためにStanとRでベイズ統計モデリングを読んでいます。

StanとRでベイズ統計モデリング (Wonderful R)

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (10件) を見る

この本の4章では，推計統計(頻度論)とベイズ統計の両方のアプローチで信頼区間と予測区間を計算してました。
これらの実装はRで行われていたため，今回はpythonを使ってこれらの計算の描写を行うことを目標にしました。*1
また，統計の定義のややこしい単語を整理しておきます。

推計統計とベイズ統計のアプローチの違い

自分の記憶の整理のために書いているので誤りがあれば指摘してください。
訂正します。

推定統計

f:id:Leo0523:20190113185610p:plain — 図1: 母集団と標本集団

推計統計では，解析対象は「母集団からランダムに抽出された標本」です。(図1)
母集団には全てのデータが存在しており，それらの平均値や分散はただ一つに決定されます。
一方で，標本は，母集団から抽出されるたびに異なるものとなるため，平均や分散は確率的なものになります。

f:id:Leo0523:20190113202346j:plain — 図2:標本集団の推定結果

図2の青色の母集団の正規分布で，黒色が標本集団から抽出した正規分布であるとします。
このように，実際の母集団の結果と，標本集団から推定した結果は微妙に異なることはありえます。
例えば，標本数Nが小さいときの推定した分散は，母集団の分散より小さくなることが計算でも証明されています。

参考 staff.aist.go.jp

N数が十分大きければ，母集団の平均と分散は推定値と完全に一致します。

ベイズ統計

ベイズ統計では，観測されたデータが全てだと視点からスタートします。
推定統計では，標本を抽出するたびに確率分布は異なりますが，ベイズ統計では図1右側の標本集団が全データです。
そのため，推定統計では母集団の平均値や分散はただ一つの値だったのですが，ベイズ統計ではこれを確率的に推定しようという考え方になります。

信頼区間と予測区間

この2つの差を途中まで理解していなかったので記載しておきます。

参考

おっと危ない：信頼区間と予測区間を混同しちゃダメ - Take a Risk：林岳彦の研究メモ

予測区間95%というのは，新しいデータは95%の確率でその区間の中に存在することを指します。*2
一方で，信頼区間は推計統計とベイズ統計で意味合いが異なります。
我々が解釈しやすいのは，ベイズ統計であり推計統計ではないことを覚えておけば

推計統計における信頼区間

推定統計において信頼区間の解釈は，「繰り返し標本を抽出したときに，その抽出した標本の信頼区間の中に母集団の平均値が入っている確率」になります。
例えば，信頼区間95%であるときには，100回標本抽出を繰り返すと95回は推定した信頼区間の中に母平均が含まれます。

下のURLの画像がわかりやすかったです。

参考 bellcurve.jp

　ベイズ統計における信頼区間

これは単純に，母集団の平均値が信頼区間に入っている確率です。
解釈がしやすいですね。

計算と可視化を行ってみた。

推定統計

github.com

コードはこちらになります。

計算はstatsmodelに行わせました。
線形回帰を行うところまでは知っていたのですが，信頼区間と予測区間は，

result = sm.OLS(Y, X).fit()
result95 = result.get_prediction().summary_frame(alpha=0.05)

としてpandasのdataframe書き出すことができます。はじめwls_prediction_stdを使っていたのですが，これは予測区間しか出すことができませんでした。

また，区間を描写するのに，matplotlibのfill_between関数を使いました。大変便利です。

f:id:Leo0523:20190113205540p:plain — 図3: 推定統計による信頼区間と予測区間。区間が細くなっている方が信頼区間

f:id:Leo0523:20190113205544p:plain — 図3: 推定統計による信頼区間と予測区間。区間が細くなっている方が信頼区間

結果として，得られたのが図3になります。これは本のものと対応できていることが確認できました。

ベイズ推定

github.com

stanの部分に関しては，

StanとRでベイズ統計モデリング (Wonderful R)

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (10件) を見る

にかかれています。pystanのAPIに関して言えば，stanファイルに記載する部分をpythonの文字列としてダイレクトに書き込み，

stanmodel='''
data {
    int N;
    real X[N];
    real Y[N];
}

parameters {
    real a;
    real b;
    real<lower=0> sigma;
} 

model {
    for (n in 1:N){
    Y[n] ~ normal(a + b*X[n], sigma);
    }
}
'''

コンパイルすることで，

sm = pystan.StanModel(model_code=stanmodel)

実行できます。今回は初めてということで，generated quantitiesを使わずに，python側で正規分布を生成しています。

結果として，

f:id:Leo0523:20190113211902p:plain — 図4:ベイズ信頼区間と予測区間

f:id:Leo0523:20190113211909p:plain — 図4:ベイズ信頼区間と予測区間

図4が得られ，無事に区間推定を行うことができました。

まとめ

pythonを使ってStanとRでベイズ統計モデリングの４章の計算を行ってみました。
信頼区間と予測区間の定義をまとめました。

*1:意外なことに観測範囲内ではあんまり誰もやってませんでした。

*2:いろいろ探ったのですが，推計統計もベイズ統計も同じように思えました。誤っていれば教えてください。

2018-12-24

reについて

python

正規表現についていつもわからなくなるので調べた関数について記録ついでに一つ一つ書きます。

正規表現の使い方を今の所よく把握していないのですが，それは追々身につけるとして，今は関数だけおっていきます。

re.sub()

note.nkmk.me

re.sub()は，第一引数に正規表現パターン，第二引数に置換先文字列，第三引数に処理対象の文字列を書きます。

これだけでだいたいわかりますね。

Python公式ドキュメントを見ると引数は5つありました。

6.2. re — 正規表現操作 — Python 3.6.5 ドキュメント

省略可能な引数 count は、置換されるパターンの出現回数の最大値です； count は非負の整数でなければなりません。もし省略されるかゼロであれば、出現したものがすべて置換されます。パターンのマッチが空であれば、以前のマッチと隣合わせでない時だけ置換されますので、 sub('x*', '-', 'abc') は '-a-b-c-' を返します。

flagsについてはよくわからなかったので，ぐぐると他の大文字，小文字の区別などをしてくれる引数のようでした。

d.hatena.ne.jp

うまい寿司が食いたい。

うまい寿司が遠慮なく食べれるようになるまで，進捗とか垂れ流すブログ

pythonで標準出力を取得する。

pythonの相対インポートがわからなくなったメモ

pytorchでpretrained modelの最終層だけlrを変更する。

PILで取り込んだ画像をnumpyに変換してmaskと重ねる

ベイズ線形回帰

動機

線形回帰の式

１次関数での線形回帰

事前分布

尤度

事前分布と尤度の積

実装

まとめ

StanとRでベイズ統計モデリング,　４章の計算をしてみた

概要

動機

推計統計とベイズ統計のアプローチの違い

推定統計

ベイズ統計

信頼区間と予測区間

推計統計における信頼区間

ベイズ統計における信頼区間

計算と可視化を行ってみた。

推定統計

ベイズ推定

まとめ

reについて

re.sub()