例えば, 以前触れた人の歩幅と距離の例でいえば, 歩幅が平均1mに対して一歩一歩が±どの程度の誤差が起きるかというのもパラメータといえそうです. ある人は±30cmでずれるかもしれませんし, ある人は±5cmかもしれません. この誤差のパラメータも結果的に何mの距離を歩けたかという観測データを見れば推定できそうです.

1.2 統計モデリングの目的

統計モデリングの目的は主に2つあります.

解釈: 現象をの理由や仕組みがわかること, 解釈できること
予測: 過去のデータから未来の振る舞いを予測すること

また, 機械学習と古典的な統計分析と比較して場合の統計モデリングの特徴です.

機械学習
- 解釈: △/✕
- 予測: ◯
古典的な統計分析
- 解釈: ◯
- 予測: ✕
統計モデリング
- 解釈: ◯
- 予測: ◯

とりあえず, 「統計モデリング最強」みたいな感じです(｀・ω・´)

1.3 確率的プラグラミング言語

確率モデルのパラメータ推定の大変さとして主に以下の2つがありました.

難解な数式を実装するのでバグりやすい
モデルが変わるたびにコードを修正しないといけない

これらの問題を解決するのが 確率プログラミング言語 で, Stanもその言語のうちの一つです.

1.4 なぜStanなのか？

確率プログラミング言語の主流としては, 以下の2つがありました.

WinBUGS
- WIndowsオンリー
- エラーメッセージがわかりづらい
- 開発が2007年から停滞
JAGS - Just Another Gibbs Sampler
- 更新されていない
- ドキュメントが少ない

そこでStanの登場.

Stan
- 2012年ころからgithubで開発開始
- R_Python_Matlabのインターフェースがある
- エラーメッセージがわかりやすい
- ドキュメントが充実
- さまざまなパラメータ推定のアルゴリズムが使える
  - NUTS: MCMCというサンプリングアルゴリズムの一種のHamiltonian Monte Carloの一実装
  - ADVI: 変分ベイズ法の一実装. NUTSより50倍速いケースがあるらしい

1.5 なぜRStanなのか？

Rは可視化やデータ加工が得意なのでいろいろ捗るらしい(｀・ω・´)

参考資料

2018-11-02

統計モデリング - pythonで『StanとRでベイズ統計モデリング』

pymc4 python r ベイズ統計モデリング stan 確率

f:id:yukinagae:20181031083449p:plain

はじめに

この本はStan（スタン）というソフトウェアとそのR用のパッケージであるRStan（アールスタン）を使って, 統計モデリングを習得する本である

イントロで書籍の目的や内容や流れ, 前提知識や対象読者が記載されていたいりするので意外と大事です. まず上記の冒頭の「何について書かれた本か？」を読んでいきます.

以下の単語について明確にしていきます. なんとなくわかったつもりになっていることも多いので, できるだけ see: [参考資料のリンク] で参考資料も調べて記載するようにします.

Stan
R
RStan
統計モデリング

統計モデリングとは

それでは以上のようなソフトウェアを用いて実施する 統計モデリング とはなんでしょうか.

統計モデリングとは確率分布を使った数理モデルをデータにあてはめて現象の理解と予測を促す営みだ.

確率分布 を使って 数理モデル をデータにあてはめて現象の理解と予測をすることということなので, 次の単語を整理してみます.

確率分布
数理モデル

確率分布

確率分布（かくりつぶんぷ, 英: probability distribution）は、確率変数の各々の値に対して、その起こりやすさを記述するものである。

see: 確率分布 - Wikipedia

Wikipediaの例によると, 例えばサイコロ2個を振ったときの出目の和は 確率変数 で, 以下のような確率の分布を 確率分布 と呼びます.

f:id:yukinagae:20181102120153p:plain

サイコロ1個の出目の範囲は1から6のため, サイコロ2個を振ったときの出目の和は2から12の範囲になります. ここで重要なのは, サイコロ2個の出目の和という確率変数があった際に, 2から12のうちどの値になりやすいかは異なるという点です.

上記の例では, （単に組み合わせが多いため）7の値になりやすいことがわかります. 一方, 2もしくは12は最も確率的に低い値となります. しかし, 思考実験として両方のサイコロが歪んでいて小さい目が出やすい場合には, 2から6の値になりやすく, 一方8〜12の値にはなりづらくなるため確率分布が変化します.

このようにある値や事柄が起きやすいかどうかの確率を確率分布と呼びます.

f:id:yukinagae:20181102120017j:plain

補足ですが, 確率については以下の本がわかりやすいです.

see: プログラミングのための確率統計

数理モデル

そもそも 数理モデル という言葉が聞き慣れないため, やっぱりWikipediaを見てみます.

数理モデルは、特に数学によって記述されたモデルのことである。モデルという言葉に含意されているように、対象とのズレ（特に近似や抽象化）が意識されていることが多い。モデルの正当性が実験や観察などによって裏付けられ、非常にうまく行っている事が確かめられている場合は「理論」と呼ばれるようになることもある。

簡単な例.

「A君が歩けば歩くほど前に進む。歩幅が広いほど前に進む。」という現象を
（距離）＝（歩幅）×（歩数）
という数式で表せば、これは数理モデルである。

see: 数理モデル - Wikipedi[f:id:yukinagae:20181102120153p:plain]a

とりあえず y = ax + b のように何らかの数式で事象を表現することを数理モデルと呼ぶようです.

つまり

統計モデリングとは, 確率分布 （というある値や事柄が起きやすいかという分布の情報）を使って 数理モデル （という数式を使って事象を説明するモデリング）をデータにあてはめて現象の理解と予測をすることと説明できます.

つまり, 先程の数理モデルの例では, （距離）＝（歩幅）×（歩数） という概念的な数式がありましたが, これに 確率分布 を取り入れると統計モデリングとなりそうです.

例えば, 人間が歩く場合常に一定の歩幅で歩くことは難しいため, ある程度の±の誤差が生まれますこれを歩幅の大きさが確率的に変動すると捉えると, 結果的に距離も確率的に変動しますと捉えることができそうです.

数理的なモデル: （距離）＝（歩幅）×（歩数）
確率的なモデル: （距離がどれくらいになりそうか）＝（歩幅の大きさがどれくらいになりそうか）×（歩数）

f:id:yukinagae:20181102120206j:plain

「じゃあ具体的にどういうこと？」という点がまだ気になってきますが, まずは概要を理解したところで次を読み進めます.

参考資料

2018-11-01

Stan/R - pythonで『StanとRでベイズ統計モデリング』

python pymc4 r stan 確率統計モデリングベイズ

f:id:yukinagae:20181031083449p:plain

はじめに

この本はStan（スタン）というソフトウェアとそのR用のパッケージであるRStan（アールスタン）を使って, 統計モデリングを習得する本である

Stan
R
RStan
統計モデリング

Stan / R / RStan

Stan - Stan は, 統計モデリングのためのプログラミング言語・ソフトウェアです（統計モデリングについては後述）

see: Stan - Stan
see: Stan - About

R: The R Project for Statistical Computing という検索しづらい名前のプログラミング言語・ソフトウェアは, 主に統計や分析で用いられます. pythonでも同様なことが可能なことが多いですが, 手軽に可視化したり, 統計値を簡単に算出したり, 分析の試行錯誤する場合にはRが適役のように思います.

see: R: The R Project for Statistical Computing
see: Rによるやさしい統計学

Rユーザはざっくり100万人はいるらしいですが, そのようなユーザに使ってもらうためにStanは RStan というインターフェースが用意されており, Rから容易にStanのコードが実行できるようです.

see: How many people use R? - Quora

またRだけでなく次のように様々なプログラミング言語のインターフェースが提供されています. 新しいところでは The Julia Language というプログラミング言語にも対応しているようです.

PyStan（Python）
MatlabStan（MATLAB）
Stan.jl（Julia）
ScalaStan（Scala）

see: Stan - Interfaces

参考資料

2018-10-31

『StanとRでベイズ統計モデリング』をpythonでやってみる

pymc4 python stan r ベイズ統計モデリング

f:id:yukinagae:20181031083449p:plain

TL;DR

今後は以下の書籍を読んで勉強した内容をブログに書いていこうと思います(｀・ω・´)

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

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

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

選定理由

現在興味を持っている ベイズ統計モデリング の書籍では, 最も網羅的かつ実践的な内容に思えるため（他の書籍だとモデルやサンプリングの結果の評価が薄い点が気になる. 実際の業務で使う際には評価の観点がとても重要になると考えるため）
pythonの本だとコードをコピペして理解した気になってしまうため, Stan/RStanという慣れていない言語の方が都合がよい
可視化や図が多くてわかりやすそう
表紙がかわいい(｀・ω・´)

やり方

30分で読む・ブログに書くが完結できる分量にする. もしそれ以上かかるようなら読む分量を分割する
RStanのコードはできるだけpythonで置き換えて実行する. コードをコピペするだけだと理解できてないことが多いが, 別言語に置き換えるにはちゃんと理解してないとできない.
当たり前だが, 自分の考えと書籍や参考資料に書いてある内容を明確に分ける

参考文献

2018-10-30

pymc4の現状整理

pymc4 python tensorflow probability

ちょっとブログ更新の間があいてたのですが、そろそろ再開したいと思います.

pymc4のソースコードリーディングがだいたい完了したので、現状のライブラリの現状を整理します.

現状整理

主な処理はtensorflow probabilityを内部的に使用
可視化はarvizを使用
他ライブラリを使いやすくするための薄いラッパー（2018/10/30現在: 56コミットしかない）
メインコントリビューターがarviz等の整備に注力しているため、pymc4自体の開発は進んでない（逆に言えば他ライブラリの整備の優先度が高い）

今後やっていきたいこと

[ ] pymc4にコントリビュートしていく（ほとんどの処理がtensorflow probabilityにあるため、主にREADMEやexampleの整備になる想定）
[ ] tensorflow probabilityのソースコードリーディング（深く理解していくために必要）
[ ] StanとRでベイズ統計モデリング (Wonderful R) を読み進めて、pymc4で実行してみる

参考資料

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

2018-10-15

サンプリングのテスト追加と微修正[6717b3e, 007cd06, d896372] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

主にサンプリングのテスト追加

コミット

2018/07/06のコミットです.

以下が変更対象ファイルです.

pymc4/model/base.py
- 不要な変数を削除
tests/test_sampling.py
- テスト追加
- pycodestyle対応

tests/test_sampling.py

正規分布の確率分布からのサンプリングをテストし, 期待通り平均0/標準偏差1になっていることを確認している.

def test_sample():
    model = pm.Model()
    
    @model.define
    def sample(cfg):
        mu = ed.Normal(0., 1., name="mu")
    trace = pm.sample(model) # => Acceptance rate: 0.9916
    assert 0. == pytest.approx(trace["mu"].mean(), 1) # => 0.014200708
    assert 1. == pytest.approx(trace["mu"].std(), 1) # => 0.9805069

pymc4/inference/sampling/sample.py

ついでに sample 関数の実装を確認する. 主に次の3つの処理を実行している.

サンプリング時の初期値を0.5で設定（乱数ではない）
tfp.mcmc （tensorflow_probabilityのMCMC実装）でサンプリング処理を定義
tf.Session 内でサンプリング処理を実行
1. サンプリング後の状態を取得
2. acceptされたサンプル数から%を計算

def sample(model,  # pylint: disable-msg=too-many-arguments
           num_results=5000,
           num_burnin_steps=3000,
           step_size=.4,
           num_leapfrog_steps=3,
           numpy=True):
    initial_state = []
    for name, (_, shape, _) in model.unobserved.items():
        initial_state.append(.5 * tf.ones(shape, name="init_{}".format(name)))

    states, kernel_results = tfp.mcmc.sample_chain(
        num_results=num_results,
        num_burnin_steps=num_burnin_steps,
        current_state=initial_state,
        kernel=tfp.mcmc.HamiltonianMonteCarlo(
            target_log_prob_fn=model.target_log_prob_fn(),
            step_size=step_size,
            num_leapfrog_steps=num_leapfrog_steps))

    if numpy:
        with tf.Session() as sess:
            states, is_accepted_ = sess.run([states, kernel_results.is_accepted])
            accepted = np.sum(is_accepted_)
            print("Acceptance rate: {}".format(accepted / num_results))
    return dict(zip(model.unobserved.keys(), states))

各々の処理を見ていく.

サンプリング時の初期値を0.5で設定（乱数ではない）

unobserved の数分の0.5の行列のリストを生成している.

    initial_state = []
    for name, (_, shape, _) in model.unobserved.items():
        initial_state.append(.5 * tf.ones(shape, name="init_{}".format(name)))

tfp.mcmc （tensorflow_probabilityのMCMC実装）でサンプリング処理を定義

    states, kernel_results = tfp.mcmc.sample_chain(
        num_results=num_results, # => 結果として必要なサンプリング数
        num_burnin_steps=num_burnin_steps, # => 最初の方の削除対象の不要なサンプリング数
        current_state=initial_state, # => 初期値（すべて0.5の行列）
        kernel=tfp.mcmc.HamiltonianMonteCarlo( # => サンプリングにはハミルトンモンテカルロ法を指定
            target_log_prob_fn=model.target_log_prob_fn(),
            step_size=step_size,
            num_leapfrog_steps=num_leapfrog_steps))

tf.Session 内でサンプリング処理を実行
1. サンプリング後の状態を取得
2. acceptされたサンプル数から%を計算

この辺りの細かい挙動は tfp.mcmc 内の実装を実際に読んでいく必要がありそう.

    if numpy:
        with tf.Session() as sess:
            states, is_accepted_ = sess.run([states, kernel_results.is_accepted])
            accepted = np.sum(is_accepted_)
            print("Acceptance rate: {}".format(accepted / num_results))
    return dict(zip(model.unobserved.keys(), states))

TODO

[ ] tfp.mcmc 辺りの実装を確認する

参考資料

Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub