Easy to type

個人的な勉強の記録です。データ分析、可視化などをメイントピックとしています。

ScipyでBrunner-Munzel検定

概要

Brunner-Munzel検定は、不当分散のときに使えるノンパラメトリック検定です。ランクから、得られたデータの中央値に有意差があるか検出します。 ノンパラメトリックの検定では、Wilcoxon-Mann-WhitneyのU検定が非常に有名ですが、この検定は不当分散の場合には第一種の過誤を得る確率が5%にならない場合があり、適切な手法として使えません。 詳しい調査には、以下のような先行記事があります。

hoxo-m.hatenablog.com

Brunner-Munzel検定

二群の平均値(代表値)の差を検定するとき

しかしBrunner-Munzel検定は、有用性のわりにどマイナーな検定なので、Rではlawstatというパッケージで利用が出来ましたが、Pythonにはパブリックで使えるいい感じのライブラリがありませんでした。 今回アップデートされたScipy Version 1.2.0でBrunner Munzel検定が使えるようになりました。使い方の記録です。

使い方

ドキュメントはこちらです。上述した奥村先生の記事と全く同じやり方にしてみます。

from scipy.stats import brunnermunzel

a = [1,2,1,1,1,1,1,1,1,1,2,4,1,1]
b = [3,3,4,3,1,2,3,1,1,5,4]
brunnermunzel(a, b)

結果

BrunnerMunzelResult(statistic=3.1374674823029505, pvalue=0.005786208666151538)

統計値やP-valueも一致しました。デフォルトでは、統計値をt分布に当てはめて推定しています。 原著論文中ではサンプルサイズが50未満であるときにはこちらを利用し、それ以上だと正規分布に対して当てはめています。正規分布に対して使いたい場合はdistributionのオプションをnormalに変更すれば使えます。 検定を両側、片側でやる場合にはalternativeが使えます。

皆様も、ぜひBM検定を使って「まだ不等分散で消耗してんの??」とドヤ顔してください。

本文

さて、この記事の目的はBM検定の宣伝でもあるのですが、真の目的はScipyにCommitする感想です。このBM検定、及びScipyのAPIドキュメント*1は僕が書きました。当該PRはこちらです。 自身のメイン言語がPythonだったこともあり、BM検定は有用そうでしたから、車輪の再開発で消耗するような人を減らしてみたかったという想いで取り組みました。 僕はOSSにコミットするのは、ほぼ初めてだったのでとてもいい勉強になりました。以下箇条書きです。

  • コード規約について
    • 基本的に、WMW検定の書き方へ似せた。ドキュメントも同じ。似た構造の検定なので、悪い戦略ではなかっただろう。
    • Scipyではgitのコミットメッセージにも規約がある。自身はそれまで非常にテキトーなコミットメッセージを書いていたので、一つの目安として自分のその他のプロジェクトでも使うようになった。
    • 自動シンタックスチェッカーを入れてFlake8で静的テストするようになった。
  • コードレビュー
    • しょーもない変数の使い方などについてもコメントをくれた。
    • 科学計算をする上で桁落ちや型変換によるエラーは大きな問題だが、それについて指摘してくれ、気がつくことが出来た。以降わりと気を使っている。
  • マージまで
    • 非常に時間がかかった。
    • 最初にコミットしたのが2017年4月だったわけだが、放置され気味でmasterにマージされるのに1年放置。
    • 何故か当初予定していたMilestone(1.1.0)から一つ先送りにされて、1.2.0でようやくマージ。結局Stagingされるまでは1年と8ヶ月かかった。
    • 1Milestoneが約半年で締められ300commit進むこと、オーガナイザーの苦労を考えると仕方ないが、もうちょっとなんとかならないのかなと思わなくはない(ガンガン貢献してオーガナイザーとして働け、というデカイまさかりが飛んできそう)。

また機会があれば、Scipyに限らずOSSに取り組んでみたいです。

*1:今読むと文法ミスが多くてしんどい。

Stanで非閉型な逆関数を含む分布のモデリング

概要

この記事はStanアドベントカレンダー-5日目の記事です(冗談です)。 趣味でやっている円周分布の勉強もそこそこ進みました。この記事では、清水本に載っている逆Batchelet変換(Inverse batschelet transformation)を使った分布をStanで実装して、現時点での問題点について考察します。

背景

円周分布として有名なものとして、von Mises分布、心臓分布、巻き込みコーシー分布などがあります。しかしこれはいずれも平均パラメータで最頻値(=モード)を取り、その左右に集中度パラメータに応じた分散が対称に存在する分布です。しかし実際のデータでは、どちらかの方向に分散が偏ることがあります。このような非対称分布の拡張は円周分布のみならず、正規分布などでも考察されています。日本語資料でしたら豊田黄色本などに記述があります。

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

さて円周分布では、次のような拡張があります。

Sine-skewed distirbution

清水本の中では、非対称な円周分布を構築する方法の一つとして正弦確率摂動法によるSine-skewed分布が載せられています。原著論文は(Abe, 2011)です。

link.springer.com

この分布は確率密度関数の設計が簡単なので、推定が容易というメリットがあります。更に確率密度関数の全区間積分値である \int_{-\pi}^{\pi}p(\theta)d\theta = 1という式が保持され、正規化が不要です。その一方で分散が集中度パラメータ(concentration parameter)だけでは制御されず、歪度パラメータからも影響を受けてしまう欠点があります。

Batschelet transformation

一方で、もともと円周分布を歪ませる方法としては(Batschelet, 1981)で提案された方法が有名です。

この方法では、 vonmises(\theta|\mu, \kappa)を基に、 aevonmises(\theta | \mu, \kappa, \nu) = vonmises(\theta + \nu\cos(\theta - \mu) | \mu, \kappa)と変換を加えます。この方法は分かりやすいし簡単なのですが、 \int_{-\pi}^{\pi}aevonmises(\theta|\mu, \kappa, \nu) d\theta = 1が維持されません。そのためvon Mises分布とは確率密度関数の形がやや変わり、 aevonmises(\theta|\mu, \kappa, \nu) = \frac{\exp(\kappa(\cos(\theta - \mu + \nu\cos(\theta - \mu)))}{c_{\mu, \kappa, \nu}}となります。 c_{\mu, \kappa, \nu}は正規化定数で、 c_{\mu, \kappa, \nu} = \int_{-\pi}^{\pi} \exp(\kappa(\cos(\theta - \mu + \nu\cos(\theta - \mu)))により求められます。確率密度関数の和が1となるように、全区間積分して和を取り、割るという至極単純なアイディアです。

しかしこの計算は積分が入って面倒です。この記事では扱いませんが、現時点ではStanに積分機能が無いので、合成シンプソン法などの数値計算で代替するのが一つのアイディアになります。なお(Batschelet, 1981)は絶版になり手に入りませんが、(Abe, 2013)でより一般的に拡張されて説明されてますので、そちらを参考にすると良いかと思います。

link.springer.com

Inverse Batschelet transformation

(Jones&Pewsey, 2012)では、逆変換 s(\theta) = t^{-1}_{1, \nu}(\theta)をかませることで、正規化が不要な非対称分布を作りました。 t^{-1}_{1, \nu}(\theta) t_{1, \nu}(\theta) = \theta -\mu - \nu(1 + \cos(\theta-\mu)逆関数です。s(t)は閉型(closed-form)で、つまり簡単な形で表すことが出来ません。このような逆関数の変換を含んだ分布の例を殆ど見ないので、これを実装してみます。参考とした論文は(Jones&Pewsey, 2012)と、Pewseyらの書いた"Circular statistics in R"です。

Circular Statistics in R

Circular Statistics in R

結果

例によってpythonです。Rの本を見てるのに…。 

逆変換で知りたい値は既知の t, \mu, \kappa, \nuのもとで t = \theta -\mu - \nu(1+\cos(\theta - \mu))を満たすような \thetaです。観測値= tが既知ということがポイントです。この場合には g(\theta) = -t + \theta -\mu - \nu(1+\cos(\theta - \mu))という関数を考えて、 g(\theta) = 0となるような \thetaを求めます。どちらにせよ非線形な問題ですが、0となる値を求めるには効率的なアルゴリズムがあるので、それを利用することができます(Wikipediaにもそんな記載があります)。3つの手法を紹介します。

シミュレーションデータは次のように作成しました。歪度パラメータが0である場合、通常のvon Mises分布と一致しますので、きちんと0が推定されるか調べます。

I = 1000
RADIAN = vonmises.rvs(size=I, loc=0, kappa=1.0)
stan_data = {
    "I": I,
    "RADIAN": RADIAN
}

Stanの組み込み関数を利用する

StanのVersion 2.17.0以降には、algebra_solverという関数で、非線形方程式の場合にも解を得ることができます。この関数はPowell法を使っています(Powell, 1970)。Stanコードの実装は次のとおりです。

  • 1~17, 48行目: 逆変換をする処理です。
  • 24~33行目: algebra_solverへ投入するデータを宣言しています。

algebra_solverはStanのマニュアルの20節に記述があります。第1引数へ解を求める関数、第2引数へその際の解、第3引数へパラメータ、第4引数へ実数型の定数値、第5引数へ整数型の定数値を使います。ちょっと使い方にはコツがあります。

  • 先程も書きましたが、今回の場合は観測値= tが既知ということがポイントです。これがパラメータになりますので、第3引数へ入れます。
  • その他モデル上でのパラメータ \mu, \nuなどもパラメータです。普通のプログラミング言語ではkargsなんかを使うのですが、第3引数にはvector型しか受け付けないので、観測値の長さより少し伸ばして代入します。
  • 定数値は今回の場合、観測値の長さしかありません。これは整数型です。実数型に該当するものは無いのですが、宣言は必要なので長さ0のベクターをtransformed_dataのブロックで宣言しておきます。

性能評価してみましょう。サンプルサイズを変えてデータを何回か生成し、真の分布に対するKL divergenceを調べます。

result = []
for I in [100, 300, 500, 1000]:
    for i in range(5):
        RADIAN = vonmises.rvs(size=I, loc=0, kappa=1.0)
        stan_data = {
            "I": I,
            "RADIAN": RADIAN
        }

        fit = model.optimizing(stan_data, init_alpha=1e-10)
        p = invBvonmises_pdf(theta, kappa=fit["kappa"], nu=fit["nu"], loc=fit["mu"])
        pt = invBvonmises_pdf(theta, kappa=1.0, nu=0, loc=0)
        kld = entropy(p, pt)
        tmp = {"I": I, "i": i, "kld": kld}
        result.append(tmp)
df = pd.DataFrame(result)

fig = plt.figure()
ax = fig.add_subplot(111)
sns.boxplot(data=df, x="I", y="kld", ax=ax)
ax.set_xlabel("Sample size")

f:id:ajhjhaf:20181126092909p:plain

データ量に応じてきちんと真の分布との距離が小さくなっていることがわかります。

ニュートン法

組み込み関数は目的関数だけあれば使える便利な手法ですが、めちゃめちゃ遅いです。ベクトル形式で渡していますが、データ数がある程度あるとほとんど使い物になりません。I=500の時点で、平均11秒ほどかかります。今回の実験がMCMCではなく、MAP推定で実行しているのもそれが理由となります。幸いにも、今回の変換は微分可能な形式で示されていますので、傾きの情報を利用したニュートン法が一つのアイディアとして考えられます。この場合には解が二次収束しますので、非常に早く計算可能です。Stanのコードは次のとおりです。

  • 1~23, 38行目: ニュートン法での逆変換です。
    • 8行目: 初期値として、 \thetaの値を使います。色々試したが、この関数の場合はこれが安定的でした。
    • 9行目:  g(\theta)を評価します。これがすべて0なら理想形です。
    • 10行目: 絶対値を取って、どれだけ0に近付ているか評価します。
    • 11行目: 反復回数の上限値を決めます。これを決めておかないと、何かあって0に近づかなかった場合は無限に計算されてしまいます。
    • 12行目: どれだけ0に近づいたら終了するか決めます。Stanの場合はmachine_precision()が計算機イプシロンとして用意されていますので、これを使います。
    • 13~15行目: tを傾きを使って更新して、誤差を再評価します。12行目で決めた基準より小さくなったら終了です。
    • 16~19行目: 反復回数が上限に達したら、反復を強制的に終了します。

先ほどと同様の誤差計算は次のとおりになります。こちらでも、問題なくKL divergenceが0に近づいています。

f:id:ajhjhaf:20181126094621p:plain

肝心の速度についてもI=500の段階で0.1秒と100倍ほど高速化することができました。2つの手法を比較すると、次のとおりになります。Powellが指数増加していることがわかります。

f:id:ajhjhaf:20181126094802p:plain

二分法

ニュートン法でハッピーエンドを迎えられるかと思いきや、そうはなりません。RstanにはStanの関数を呼び出すexpose機能があります。 s(t)がどのような挙動をするか調べてみましょう。先程のモデルを適当に保存してください。

library("rstan")
rstan::expose_stan_functions("Desktop/invBvonmises.stan")
x = seq(-pi, pi, length=1000)
mu = 0
nu = 0.99
I = 1000
y = inv_transformation_Batschelet(x, mu, nu, I)
plot(x, y)

f:id:ajhjhaf:20181126095608p:plain

良さそうに見えます。が、ここでnuの値を極めて大きくしてみましょう。

f:id:ajhjhaf:20181126095532p:plain

ぶっ壊れた値が出てきてしまいました。これは、 \nuが大きい場合に、勾配が消失することが原因と思われます。 g(\theta) \nuが大きい場合に書いてみますと、次のとおりです。

mu = 0
nu = 0.99
theta = 0
x = np.linspace(-np.pi, np.pi, 1000)
y = x - nu * (1 + np.cos(x - mu)) - theta

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y)
ax.set_xlabel("theta")
ax.set_ylabel("g(theta)")

f:id:ajhjhaf:20181126095943p:plain

thetaが-2から-1の場合にかけて、勾配が消失していることが見て取れます。この場合には、傾きを使わない手法として、二分法を使うことが考えられます。Stanのモデルは次のとおりです。

  • 1~34行目: 二分法による逆変換です。
    • 10-11行目: 上界、下界を決めます。求める解の範囲が -\piから \piなので、それぞれ g(\theta)が0以上、以下となるように決めましょう。
    • 12行目: 上界、下界から真ん中の値をとります。
    • 13-14行目: ニュートン法と同様に、 g(\theta)の評価をして、誤差を調べます。
    • 17-23行目: 評価値の正負に応じて、上界、下界を狭めます。pythonだったらnumpyのsmart indexを使えば一行でかけますが、Stanだとfor文で一つ一つ判定しなくてはいけません
    • 28行目: 二分法は一次収束です。ニュートン法と同等の精度を得るために、1000(>302)回を反復上限とします。

R側から呼び出して、同様に計算させてみると、変換できたことがわかります。

f:id:ajhjhaf:20181126104135p:plain

が、ここで大きな問題があります。R側での動作を見るに動きそうなこのモデルですが、実際に動かしてみるとLS failedが発生してうまくいきません。ナンデ……。大きな \nuから作ったデータに対してフィッティングしてエラーが起こるなら、変換の勾配が消失して微分できず失敗すると納得できるのですが、ここまで使ってきた \nu=0のデータに対して実行してもやはり失敗します。

まとめ

陽的に式をかけない変換を含むモデリングを実行することができました。Powell法だと(おそらく)数値的に安定なので使いやすいですが、データが多い場合には変換が非効率でめちゃめちゃ時間がかかります。ニュートン法は超絶早いですが、傾きがある場合じゃないと使えません。手元でまともに動く範囲を調べて、パラメータ制約つけるのがいいでしょうか…。二分法はなぜかうまくいきません。困った(ちなみにdiscourseで絶賛助けを求めています、コメント歓迎です……)。

今回利用したコードは、こちらにあります。再現実験したい方はどうぞ。

謝辞

このブログ記事を書くにあたり、StanのdiscourseでBen Bales氏に助けていただきました。お礼を申し上げます。

StanのThreading動作確認

TL; DR

導入

PyStanは2.18.0系からStanに搭載されたThreading機能を使うことができるようになりました。公式のドキュメントはこちらになります*1

この機能は、今までのMCMC chainごとにスレッドをフォークして並列化するのではなく、MCMCサンプル内での並列化です。先程のリンク先のExampleの欄にも書かれているとおりに、chain数×設定スレッド数分の並列化ができるように成りました。ただし、まだかなり作りかけの機能であるため、運用上は幾つか疑問点があります。

  • 並列化することでどれほど早くなるのか?
  • 単一のMCMC chainの場合でも並列化されるのか?
  • 並列化したことで推論結果が悪くなることはあるのか?
  • 設定スレッドはモデルのコンパイル時のみ影響するのか?サンプリング時にも影響するのか?

以上の疑問を解決するために、モデルを実際に動かしてみました。

手法

環境

次のとおりです。

  • pystan 2.18.0
  • python 3.6.5
  • CPU数 32
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from scipy.stats import norm
from pystan import StanModel
import pystan

print(pystan.__version__)
print(os.cpu_count())

モデル

Exampleに記載されているeight-school problemのデータとモデルをそのまま利用します。

model_code = """
functions {
  vector bl_glm(vector mu_sigma, vector beta,
                real[] x, int[] y) {
    vector[2] mu = mu_sigma[1:2];
    vector[2] sigma = mu_sigma[3:4];
    real lp = normal_lpdf(beta | mu, sigma);
    real ll = bernoulli_logit_lpmf(y | beta[1] + beta[2] * to_vector(x));
    return [lp + ll]';
  }
}
data {
  int<lower = 0> K;
  int<lower = 0> N;
  vector[N] x;
  int<lower = 0, upper = 1> y[N];
}
transformed data {
  int<lower = 0> J = N / K;
  real x_r[K, J];
  int<lower = 0, upper = 1> x_i[K, J];
  {
    int pos = 1;
    for (k in 1:K) {
      int end = pos + J - 1;
      x_r[k] = to_array_1d(x[pos:end]);
      x_i[k] = y[pos:end];
      pos += J;
    }
  }
}
parameters {
  vector[2] beta[K];
  vector[2] mu;
  vector<lower=0>[2] sigma;
}
model {
  mu ~ normal(0, 2);
  sigma ~ normal(0, 2);
  target += sum(map_rect(bl_glm, append_row(mu, sigma),
                         beta, x_r, x_i));
}
"""

stan_data = dict(
    K = 4,
    N = 12,
    x = [1.204, -0.573, -1.35, -1.157,
         -1.29, 0.515, 1.496, 0.918,
         0.517, 1.092, -0.485, -2.157],
    y = [1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1]
)

パラメータ

モデルは次の4つをコンパイルしました。

  • extra_args等の一切を無指定(以下Vanillaと呼称)
  • extra_argsを指定 + STAN_NUM_THREADS=1と指定(THREADS1と呼称)
  • extra_argsを指定 + STAN_NUM_THREADS=4と指定(THREADS4と呼称)
  • extra_argsを指定 + STAN_NUM_THREADS=8と指定(THREADS8と呼称)

これとは別に以下のパラメータをサンプリング時に指定しました。

  • サンプリング直前にSTAN_NUM_THREADSを再設定
  • n_jobs(公式APIでは"Sample in parallel. If -1 all CPUs are used. If 1, no parallel computing code is used at all, which is useful for debugging."と書かれているが、これがThreadingに影響するのか、MCMC chainだけに影響するのか不透明)
  • iter(MCMCの反復数)
  • chains(MCMCのchain数)

実験

上述のパラメータを設定した上で、jupyter labの%timeitオプションを付けて実行しました。各パラメータセットごとに7回反復されて、平均が取られました。notebookはこちらです。

結果

index compile num_threads sampling num_threads n_jobs iter chains avarage time std time
1 Vanilla 1 4 2000 2 0.883 0.121
2 Vanilla 1 2 2000 2 0.834 0.0624
3 Vanilla 1 4 2000 2 0.935 0.181
4 Vanilla 1 1 2000 2 1.21 0.413
5 Vanilla 1 1 2000 1 0.554 0.0567
6 Vanilla 1 1 20000 1 4.36 0.296
7 Vanilla 1 8 2000 2 0.82 0.0896
8 1 1 4 2000 2 0.936 0.0743
9 1 1 2 2000 2 0.967 0.101
10 1 1 4 2000 2 0.937 0.0693
11 1 1 1 2000 2 1.27 0.0501
12 1 1 1 2000 1 0.678 0.0833
13 1 1 1 20000 1 5.34 0.804
14 1 1 8 2000 2 0.962 0.0777
15 4 4 4 2000 2 3.63 0.547
16 4 4 2 2000 2 3.54 0.434
17 4 1 4 2000 2 0.99 0.0883
18 4 4 1 2000 2 5.4 0.576
19 4 4 1 2000 1 2.68 0.216
20 4 4 1 20000 1 21.1 1.11
21 4 4 8 2000 2 3.27 0.274
22 8 8 4 2000 2 3.28 0.456
23 8 8 2 2000 2 3.43 0.263
24 8 1 4 2000 2 0.865 0.13
25 8 8 1 2000 2 4.91 0.434
26 8 8 1 2000 1 2.52 0.314
27 8 8 1 20000 1 22.8 4.12
28 8 8 8 2000 2 3.31 0.341

average timeとstd timeの単位は秒です。

  • VanillaはSTAN_NUM_THREADSが1の場合と殆ど変わりません。中途処理があるからか、0.1秒ほど遅くなっています(Compare 1, 2, 7... vs 8, 9, ..., 14)。

  • n_jobsはchain数より大きい値を入れても早くなりません(Compare 15 vs 16等)。

    • 一方で、chain数より小さな値を入れれば、並列でMCMCできなくなるので遅くなります(Compare 15 vs 18等)。
    • また、STAN_NUM_THREADSのパフォーマンスにも影響しないようです(Compare 22 vs 23)。
  • 環境変数に設定したSTAN_NUM_THREADSは、サンプリング時にも影響するようです(Compare 15 vs 17等)。

    • コンパイル時の値は上限スレッド数というだけで、実効値は環境変数を参照するのでしょう。
    • 即ち、SGEなどのqueuing systemを使っているときには、サンプリング時に環境変数をセットする必要があります。
  • 推論結果は、見た感じ悪化しませんでした(実行結果を参照)。

  • 一番大事なことですが、 Threadingをしても遅くなるだけです (Comapare 2 vs 9 vs 16 vs 23)。

    • モデルが単純でデータが少ないからでしょうか?iter数を増やしてみればThreadingの恩恵があるかと思いましたが、傾向は変わりませんでした(Compare 6 vs 13 vs 20 vs 27)。
    • ぜひどなたかに、複雑なモデルで実験してほしいです(自分には適当なものが思いつきませんでした)。
    • %timeitがマルチスレッドの場合は時間×スレッド数をはじき出すかと思い、timeモジュール+自分の時間間隔でも計測しましたが、やはり遅かったです。
  • サンプリング時にhtopコマンドでCPU利用数を見ていましたが、高々200%ぐらいになるだけで、本来理想としていた800%等になることはありませんでした。


追記(2018-10-10)

@nan_makersstat さんや id:StatModeling (@hankagosa)さんから次のような指摘をいただきました。

というわけで、軽く実験してみます。モデルは簡単な一次回帰です。シミュレーション的にデータを10万点作成します。 注: 以下のモデルはThreadingで必要なmap_rect関数が使われていないので早くなるはずはありません。map_rect用のモデルは後述しました(2018-10-19)

model_code = """
  data {
    int I ;
    real x[I] ;
    real y[I] ;
  }
  
  parameters {
    real a ;
    real b ;
  }
  
  model {
    for (i in 1:I){
      y[i] ~ normal(a * x[i], b) ;
    }
  }
  
  generated quantities {
    real log_likelihood[I] ;
    for (i in 1:I){
      log_likelihood[i] = normal_lpdf(y[i] | a * x[i], b) ;
    }
  }
"""

I = 100000
a = 1
b = 1
x = np.linspace(-1, 1, I)
y = norm.rvs(loc=a*x, scale=b)

stan_data = {
    "I": I,
    "x": x,
    "y": y
}

これをVanilla, STAN_NUM_THREADS=1, STAN_NUM_THREADS=4で投入してみました。計測方法は同じで7回の平均です。

index compile num_threads sampling num_threads n_jobs iter chains avarage time std time
1 Vanilla 1 2 100000 2 111 5.5
2 1 1 2 100000 2 113 5.65
3 4 4 2 100000 2 108 5.71

うーん、ちょっとだけ早くなっています。しかしt検定で有意差が出る気がしません。おそらくデータ数を増やせばより効果があるのではないでしょうか。 何れにせよ、小標本なのにThreadingが必要なほどサンプリングで時間がかかりすぎている場合は、モデルが適切でない可能性をまず考えるべきです。丁寧に、メカニズムを考えて、モデルを改良していきましょう!


追記(2018-10-18)

discourseで質問したところ、ThreadingやMPIのような操作をするためには、データやモデルをそれ専用に作り変える必要があるとの回答を頂きました。 まだ試せていませんが、cmdstanを使った例についてはこのリポジトリが参考になります。ご参考までに。


追記(2018-10-19)

map_rect用のモデルを別に書いて実験しました。ついでにデータ数は合成に100万まで増やしてみました。実験につかったnotebookはこちらです。

index compile num_threads sampling num_threads n_jobs iter chains avarage time std time
1 Vanilla 1 1 500 1 128 6.71
2 4 4 4 500 1 173 5.71

なんと、やっぱり遅くなってしまいました。しかし、やはりtop等のコマンドで確認してもCPU使用率が100%までいかないんですよね。別の原因がある気がしますが、現状不明です。

*1:この機能は単体のマシンでの並列化であり、MPI 環境を使った並列化はまた別に実装されています

離散ハート型分布の導出とシミュレーション

概要

相変わらず、ちょこちょこと円周統計を勉強しています。

テキスト中では、方向統計の分布としてハート型分布が紹介されています。この分布は

 f(\theta) = \frac{1 + 2\rho\cos{(\theta)}}{2\pi}

確率密度関数を持つ分布です。確率密度関数や、そこからサンプリングした乱数を可視化すると次のようになります。左側の図がユークリッド空間、右側の図が極座標系です。

f:id:ajhjhaf:20180724124401p:plain

他にも円周分布として、有名なvon Mises分布、巻き込みコーシー分布、巻き込み正規分布、一般化ハート型分布(or Jones-Pewsey分布)などなども書かれています。

以上に挙げた分布は全て連続分布です。テキスト中では一節割かれて、格子分布が紹介されていますが詳しくありません。しかし実際に角度のデータが観測された場合には、有限個の離散値になることが多いのではないでしょうか。この記事ではハート型分布を離散化することで、分布を導出し、乱数生成もしてみます。

分布の導出

先程も書いたハート型分布を変形させて、離散化させてみます。アプローチとしては単純で、確率分布 f(\theta)に対して、

 \int_{\theta-\alpha}^{\theta+\alpha} f(\theta) d\theta

を計算します。円周全体が L個に分割されているとしたら、それぞれの \alphaは前後の角度との差分になるので、 \alpha = \frac{2\pi}{2L} = \frac{\pi}{L}で求まります。というわけで、以上の結果をSympyで計算してみます。

import sympy
init_printing(use_unicode=True)
from sympy import symbols, integrate, simplify, init_printing, cos, pi
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all'

theta, mu, rho, Theta, alpha = symbols("theta, mu, rho, Theta, alpha", real=True)
c_pdf = (1 + 2 * rho * cos(theta - mu)) / (2 * pi)
c_pdf_int = integrate(c_pdf, (theta, Theta-alpha, Theta + alpha))
simplify(c_pdf_int)

結果は  \frac{1}{\pi} \left(\alpha + 2 \rho \sin{\left (\alpha \right )} \cos{\left (\Theta - \mu \right )}\right) = \frac{1}{\pi} \left(\frac{\pi}{L} + 2 \rho \sin{\left (\frac{\pi}{L} \right )} \cos{\left (\Theta - \mu \right )}\right)

として得られます。とりあえず適当に、離散ハート型分布などと名付けます*1

乱数の生成

どのような挙動をするか興味がありますので、この分布に従った乱数を生成してみます。乱数生成についても上述したテキストに書いてあります。よく使われるのは逆関数法ですが、ハート型分布では逆関数の導出が難しいので現実的ではないそうです。離散ハート型分布も同じでしょう。ここでは並べて紹介されている棄却採択法(棄却サンプリング)を使います。といっても、テキストで紹介されていた棄却採択法は連続分布に対するものだったので、MathWorksの記事を参考にしました。手順は殆ど変わらず、乱数を作りたい分布 P_p(\Theta)に対して

  1. 任意の分布 P_q(\Theta)を選ぶ
  2.  P_p P_qからサンプリングされる任意の値 p_i,  q_iに対して、 \frac{p_i}{q_i} \leq cを満たすようなcをとる
  3. [0, 1]の一様乱数から乱数 uを生成する
  4.  P_qから乱数 vを生成する
  5.  u \leq \frac{P_p(v)}{cP_q(v)}を満たすようなら乱数として採択、満たさないなら棄却する
  6. 3に戻る

を繰り返すだけです。

テキストでは連続分布であったため、 P_qとして \frac{1}{2\pi}を使っていましたが、この問題では離散分布であるので、 \frac{1}{L}を使います。このときに

 \frac{1}{\pi} \left(\frac{\pi}{L} + 2 \rho \sin{\left (\frac{\pi}{L} \right )} \cos{\left (\Theta - \mu \right )}\right) \leq \frac{c}{L}を満たすようなcの最小値は

 c = 1 + \frac{2L}{\pi}\rho\sin{(\frac{\pi}{L})}

となります。単純に不等式を解くだけです。これを実装してみると次のようになりました。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from collections import Counter
from scipy.stats import randint

def polar_twin(ax):
    ax2 = ax.figure.add_axes(ax.get_position(),
                             projection='polar',
                             label='twin',
                             frameon=False,
                             theta_direction=ax.get_theta_direction(),
                             theta_offset=ax.get_theta_offset())
    ax2.xaxis.set_visible(False)

    for label in ax.get_yticklabels():
        ax.figure.texts.append(label)

    return ax2

def dcardioid_pmf(theta, loc, rho, L):
    pmf = (np.pi / L + 2 * rho * np.sin(np.pi / L) * np.cos(theta - loc)) / np.pi
    return pmf

# AR法での実装
def dcardioid_rvs(rho, loc, size, L):
    c = 1 + 2 * L * rho * np.sin(np.pi / L) / np.pi
    result = []
    while len(result) != size:
        u = np.random.uniform(0.0, 1.0)
        v = randint.rvs(0, L) * 2 * np.pi / L
        p = dcardioid_pmf(v, rho=rho, loc=loc, L=L)
        q = 1 / L
        if u < p / (c * q):
            result.append(v)
    return result

# 乱数の長さ
I = 5000
# 区分する数(L)
num_r = 360
#分布のパラメータ
mu = 0
rho = 0.45

# 範囲
r = np.arange(0, 2 * np.pi, 2 * np.pi / num_r)
# 確率密度関数の取得
f1 = dcardioid_pmf(r, loc=mu, rho=rho, L=num_r)
# 乱数生成
y = dcardioid_rvs(rho=rho, size=I, loc=mu, L=num_r)
# 乱数をgroupingして数え上げる
c = Counter(y)
ydf = pd.DataFrame(list(c.items()), columns = ["theta", "count"])
ydf = pd.DataFrame(np.arange(0, num_r, 1) * 2 * np.pi / num_r, columns = ["theta"]).merge(ydf, how="outer", on=["theta"])
y_hist = ydf["count"]
width = 2 * np.pi / num_r


# 可視化
fig = plt.figure(figsize=(10,5))

# ユークリッド系
ax = fig.add_subplot(1,2,1)
ax.plot(r, f1)
ax2 = ax.twinx()
ax2.plot(r, y_hist)

# 極座標系
ax = fig.add_subplot(1,2,2, projection="polar")
ax.plot(r, f1)
ax.set_theta_zero_location("N")
ax.set_rticks(np.linspace(0, round(ax.get_rmax()+0.05, 1), 3))
ax2 = polar_twin(ax)
ax2.bar(r, y_hist, width=width*0.90)
ax2.set_theta_zero_location("N")
ax2.set_rticks(np.linspace(0, round(ax2.get_rmax(), 0), 3))
ax2.set_rlabel_position(ax.get_rlabel_position() + 180)
  • 極座標系で二軸プロットをするためには、以下のQ&Aを参考に改造しました。
    • 最新版のmatplotlibにはラベルを指定するメソッドが追加されているので、polar_twinの中でラベルの角度を指定する必要がありません。

stackoverflow.com

  • 離散角度の数え上げはcounter + pandasでjoinという超まどろっこしいことをしています。もっと美しいコードはあるはずです。

    • np.histogramでやってみたのですが、境界値の選択が雑なのか、上手くいきませんでした。
  • 棄却採択法での乱数生成ももっと綺麗に書けますが、上述した解説と対応させました。

結果は次のとおりになります。シミュレーション結果も上々です。

f:id:ajhjhaf:20180724125649p:plain

他の分布でも使えるか?

今回はハート型分布を離散化することが出来ましたが、この手順による拡張は他の分布にも使うことが出来るでしょうか?個人的には結構難しいと思います。

巻き込みコーシー分布を考えます。同じように積分すると、

 \begin{cases} - \frac{1}{4 \pi} \left(\rho^{2} - 1\right) \left(\tan{\left (- \frac{\Theta}{2} + \frac{\alpha}{2} + \frac{\mu}{2} \right )} + \tan{\left (\frac{\Theta}{2} + \frac{\alpha}{2} - \frac{\mu}{2} \right )}\right) & \text{for}\: \rho = -1 \\\frac{1}{4 \pi} \left(\rho^{2} - 1\right) \left(\frac{1}{\tan{\left (\frac{\Theta}{2} + \frac{\alpha}{2} - \frac{\mu}{2} \right )}} + \frac{1}{\tan{\left (- \frac{\Theta}{2} + \frac{\alpha}{2} + \frac{\mu}{2} \right )}}\right) & \text{for}\: \rho = 1 \\\frac{i}{2 \pi} \left(\log{\left (\frac{1}{\rho + 1} \left(- i \rho + \left(\rho + 1\right) \tan{\left (- \frac{\Theta}{2} + \frac{\alpha}{2} + \frac{\mu}{2} \right )} + i\right) \right )} - \log{\left (\frac{1}{\rho + 1} \left(- i \rho - \left(\rho + 1\right) \tan{\left (\frac{\Theta}{2} + \frac{\alpha}{2} - \frac{\mu}{2} \right )} + i\right) \right )} - \log{\left (\frac{1}{\rho + 1} \left(i \rho + \left(\rho + 1\right) \tan{\left (- \frac{\Theta}{2} + \frac{\alpha}{2} + \frac{\mu}{2} \right )} - i\right) \right )} + \log{\left (\frac{1}{\rho + 1} \left(i \rho - \left(\rho + 1\right) \tan{\left (\frac{\Theta}{2} + \frac{\alpha}{2} - \frac{\mu}{2} \right )} - i\right) \right )}\right) & \text{otherwise} \end{cases}

が得られました。わお。これを使うのはしんどいでしょ。

von Mises分布では

 \frac{1}{2 \pi I_{0}\left(\rho\right)} \int_{\Theta - \alpha}^{\Theta + \alpha} e^{\rho \cos{\left (\mu - \theta \right )}} d\theta

が得られます。この積分は陽に解析解が出ないようです。

積分する方法以外に、各分布の分布関数を \alphaだけ変化させたときの差分を確率密度関数として使うことも考えましたが、巻き込みコーシーもvonMisesでも式が綺麗にならなかったので諦めました*2。連続的に便利な分布でも、簡単に離散化できないんだなぁ、といういい勉強になりました。

*1:discrete cardioid distributionで調べると一件だけ統計数理研究所のセミナーがヒットしますが、他に資料がありません。参考にしたいなぁ。。

*2:というか、von Mises分布は恒等関数や無限和が入ってきて扱いづらい

角度データの変化検知を可視化するプロット

概要

角度データの解析を勉強しています。

この書籍中ではP95で、時系列角度データの変化点の検出について1節割かれています。 しかし実際に可視化した際にどのような挙動を示すのか図表がなく分からなかったので、可視化してみました。

データの生成

いつも通りjupyterで実験しました。 簡単にvon mises分布を2つ容易して、それを結合します。以下の例ですと系列店700の段階で傾向が変わったということですね。可視化した際にもなんとなくわかりますが、von mises分布自体のノイズが大きいので判断が難しいです。

import matplotlib.pyplot as plt
from scipy.stats import vonmises
%matplotlib inline

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
I = 1000
i1 = 700
i2 = 300
x = np.arange(0, I)
y1 = vonmises.rvs(kappa=0.8, loc=0, size=i1)
y2 = vonmises.rvs(kappa=0.4, loc=np.pi/2, size=i2)
y = np.concatenate([y1, y2])
ax.plot(x, y)

f:id:ajhjhaf:20180418201352p:plain

累積和プロット(Cumulative sum plot; CUSUM plot)

CUMSUMプロットでは、生成された角度を正弦と余弦に分解した時に、それぞれの成分の累積値をプロットします。同じ分布から生成された角度の累積値は一定の傾向を持ち続けますが、分布が変化した時にこの傾向が変わる、といった点から判断できるようです実際にノイズが取り除かれて、可視化結果も分かりやすくなっています。

fc = np.vectorize(lambda x: np.cos(x / 2 / np.pi))
cosy = fc(y)
c = np.cumsum(cosy)

fs = np.vectorize(lambda x: np.sin(x / 2 / np.pi))
siny = fs(y)
s = np.cumsum(siny)

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(c, s)
ax.set_xlabel("Cumulative cosine value")
ax.set_ylabel("Cumulative sine value")

f:id:ajhjhaf:20180418201929p:plain

累積平均方向プロット(cumulative mean directional plot)

このプロットでは、系列における時点jまでの累積平均方向をプロットします。分布が変化すると、平均の方向も変化するのでやはり判定できるといったものです。テキストではcos, sinの累積成分 C_j,  S_jとそこから求められる合成ベクトルの長さ R_jに対して

 \cos{\bar{\theta_j}} = \frac{C_j}{R_j}

あるいは

 \sin{\bar{\theta_j}} = \frac{S_j}{R_j}

を満たす \theta_jをプロットすると書いてあります。ですが、これを求める際には

 \tan{\bar{\theta_j}} = \frac{\sin{\bar{\theta_j}}}{\cos{\bar{\theta_j}}} = \frac{\frac{S_j}{R}}{\frac{C_j}{R}} = \frac{S_j}{C_j}

なので、別に R_jは必要ありません。また系列の最初の方ではnが少なくて値が安定しません。なので、最初100を削ってスケールを合わせたものも出してみました。

y_bar = np.arctan(s/c)

fig = plt.figure()

ax = fig.add_subplot(2,1,1)
ax.plot(x, y_bar)
ax.set_xlabel("series")
ax.set_ylabel("average direction")

ax = fig.add_subplot(2,1,2)
ax.plot(x[100:], y_bar[100:])
ax.set_xlabel("series")
ax.set_ylabel("average direction")
ax.set_xlim(0, I)

f:id:ajhjhaf:20180418202829p:plain

ランク累積和プロット(rank CUSUM plot)

ランク化して計算します。今回扱うプロットの中では一番分かりやすいですね。元々の分布がなんだったかにも影響しそうですが、計算が面倒くさいだけあっていい感じです。

from scipy.stats import rankdata
gamma = 2 * np.pi * rankdata(y) / y.size
U = np.sqrt(2 / y.size) * np.cumsum(np.cos(gamma))
V = np.sqrt(2 / y.size) * np.cumsum(np.sin(gamma))
B = np.sqrt(np.power(U, 2) + np.power(V, 2))

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(x / I, B)
ax.set_xlabel("relative series")
ax.set_ylabel("B")

f:id:ajhjhaf:20180418203106p:plain

累積円周分散プロット(cumulative circular variance plot)

テキスト中では

各j=1,...,nに対して、円周分散1-R_jを計算し…

と書いてありますが、正しくは合成ベクトルの大きさ R_jではなく、平均合成ベクトル長 \bar{R_j}を使って計算する必要があります。傾向が累積平均方向ベクトルと変わらないですね。

r = np.sqrt(np.power(c, 2) + np.power(s, 2))
r_bar = r / np.arange(1, I+1)
v = 1 - r_bar

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(x, v)
ax.set_xlabel("series")
ax.set_ylabel("cumulative circular variance")

f:id:ajhjhaf:20180418203433p:plain

まとめ

可視化を駆使することで、時系列で得られた角度データの変化を分析することが出来ました。分析に役立てていきましょう。

Python + Pyomoによる(非線形)数値最適化

TL; DR

  • Pythonのデータマネジメント技術と数値最適化をスムーズに繋げたい
  • Pyomoを使うことで自然な記法でモデルを組み立てることが出来る
  • Webドキュメントは貧弱だが、コミュニティは活発!

概要

数値最適化は機械学習や数値モデリングの基礎も基礎ですが、単に役に立ちます。Pythonという観点を度外視すれば、いろいろな手段がありますが、基本的には共通するフローがあります。

  • 現象から最適化したい目的関数を決定する
  • 目的関数に付随する制約条件を立式化する
  • データを用意し、モデリング言語を記述する
  • モデルをソルバーへ投入し、解を得る

こういった作業をするための環境は商用でも色々と売られていまして、例えばAIMMS社とか、GAMS社AMPL社といった会社は最適化の開発環境や言語を作って販売しています。

もちろん非商用でも、こういう団体は存在し、COIN-ORではライブラリやソルバーを開発して公開しています。商用のソフトウェアに比べれば少々パフォーマンスに劣りますが無料で使えるということを考えれば文句も出ません。

しかしながら、例えばAMPL言語の書き方を見てみると…ちょっとめんどくさそうです。Stanかよ。新しく言語を覚えるのは大体の場合面倒くさいので、実務でちゃちゃっと最適化したい場合には向きません。どうせだったらPythonで済ませられたらうれしいです。

Pythonで最適化?

Pythonで最適化する方法も色々あります。

Scipy.optimize

一番最初に思いつく方法だと思います。アルゴリズムも色々ありますし、導入も簡単です。ただし、簡単な最小化とか線形計画法(LP)に特化しています。もうちょっと痒いところに手が伸びるものが欲しいのです。

PuLP

日本語で調べるとこのライブラリについての記事が一番多く出てきます。PuLPは先程のCOIN-ORによって開発が進められているライブラリです。使いやすいですし、Sphinxでのドキュメントも有るのが個人的にはありがたいところです。

ただし自分が使う上では、

  • 使えるソルバーが少ない
  • 最近開発が停滞気味

といった2点が気になりました。

リポジトリに書いてあるように、PuLPはLP/MILP用ライブラリです。デフォルトソルバーはCOIN-OR CBCで、他に使えるのもCPLEXGLPKなどLPソルバーです。非線形問題は解けません。

もう一つ、PuLPは2017年の中頃からレポジトリが更新されていません。ドキュメントも更新がされておらず、実際に使う上で良く分からないことも多いです。というか英語ドキュメントすらなさすぎて、日本語記事のほうが詳しいかも…。先行記事では

などがありました。他にもQiitaやBlog上に、沢山の詳しい記事があります。

CVXPY

こちらの記事がとても詳しいです*1

myenigma.hatenablog.com

Pyomo

PyomoはOpen-OR のウェブページにも載っていますが、メインで開発しているのは米サンディア国立研究所のようです。昔はCooprという名前でした。Pythonオブジェクト指向な書き方をとても重視したライブラリです。しかし僕はPuLPの次にこれを使ってみたら、結構似ていると感じたので、移行も簡単でしょう。PyomoやPuLPは2008年ぐらいに沢山開発されたPythonライブラリで、今も生き残っているものたちです。なので、書きやすさは担保されていると思います(参考↓)。

Pyomoの嬉しいところは他に2つあって、

  • 使えるソルバーが多い
  • コミュニティがそこそこ活発

という点です。PyomoはNLフォーマット)を利用するソルバー(AMPL Solver Library)に対応しています。PuLPで使えていたCBC、GLPK、CPLEXなども利用することが出来ますし、その他についてはAMPLのページを見たほうが早いですね。PuLPでは使えなかった非線形ソルバーで有名なSCIPが使えますし、Open-OR手動で開発されているBonminやCouenne、IPOPTなども利用できます。ありがたいものです。

もう一つ、PyomoはPuLPと同じぐらいドキュメントは弱いのですが*2、コミュニティは活発です。リポジトリは頻繁にコミットがありますし、チュートリアル用のリポジトリモデル例もあります。個人的に一番ありがたかったのは、Google Groupが活発で、質問にすぐに答えてくれた点です。

そういうわけで、今回はPyomoを使って最適化を試してみます。先行記事には

がありました。他に日本語の記事があまり出てこないので、残しておきます。

導入

ライブラリ自体は

$ pip install pyomo

で終わりですが、これだとソルバーが全然入っていません。ソルバーはコマンドラインから実行できるようにする必要があります。自分はOSXの端末に以下のソルバーを入れました。皆様はお好きにどうぞ。

  • glpk
    • brew install glpkで終了!
  • CPLEX
    • アカデミックライセンスを使いました。
    • IBMインストーラーで導入後、/Applications/CPLEX_Studio128/cplex/bin/x86-64_osxへPATHを通せばOKです。
  • IPOPT
    • こちらを参照しました。少し古い記事ですが、この通りで動きます*3
  • SCIP
  • Bonmin、Couenne
    • Open-ORのソルバーなので、IPOPTと同様の手順で入ります。
    • IPOPTの時にOPENBLASやLAPACKを入れていれば、飛ばしてOKです。

非線形最適化問題を解く

どうせなので、非線形問題を解いてみましょう。例として、放送大学の資料にあった問題を解きます。

工数の最適配分の計算

* 工程i(i=1, 2, ..., n)にx_i(日)かけたときの仕上がりの程度がlog(a_i * x_i + 1) ただしa_i > 0 で決定する
*  各工程には最低でt_i(>0)日必要である
*  全行程をT日以内に終わらせる必要がある
*  仕上がりを最大化するには、各工程にどのように割り振ればいいか

仕上がりの程度にlogが入っているので非線形最適化問題になります。細かく言うなら、変数であるxが整数しか取らないので、整数非線形計画法(Integer Non Linear Programming; INLP)というパラダイムで解かれる問題ですね。ただし整数計画問題(Integer Programming; IP)やより狭義の整数線形計画問題(Integer Linear Programming)に特化した研究は色々ありますが、INLPに特化したものは少なく、一部の変数だけが整数条件を取る混合整数非線形計画法(Mixed INLP; MINLP)の一部として取り組むことが多いです。ソルバーもMINLP用のものを探したほうがいいです。

以下Jupyter + Pythonでの解法です。

  • 2, 26行目: これはJupyter環境でやっています。Pyomoの例題ではコマンドラインから動かしていますが、こうすればJupyterでも動きます。optの部分でsolverを指定して下さい。
  • 7, 9行目: PyomoはList型のIndexingを自動でやってくれないので、Dict型に変換しています。
  • 14行目: 変数についてはVarで宣言し、2つめの位置変数へデータ型を入れます。今回は日数なので、非負整数です。
  • 16~20行目: 制約条件を書いていきます。最適化条件と順番が前後しても大丈夫です。
    • Pyomoの変な特性ですが、グローバル変数とローカル変数がごっちゃになっても受け付けてしまいます。より良い書き方がありそうですが…現状知りません。
  • 21行目: 制約条件時に、Iの各要素iで成立…としたい場合は、第一位置変数へ入力すれば受け付けてくれます。
  • 22~24行目: 最適化条件を書きます。senceで最大化(maximize) or 最小化(minimize)を指定します。
  • 27行目: opt.solveをする時に、tee=Trueとすればソルバーのメッセージが表示されます。デバッグ時に便利です。
  • 28行目: モデルの結果を表示します。

以下出力です(デバッグ部分は除きました)。

Model Days allocation

  Variables:
    x : Size=4, Index=I
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  15.0 :  None : False : False : NonNegativeIntegers
          2 :     0 :  20.0 :  None : False : False : NonNegativeIntegers
          3 :     0 :  50.0 :  None : False : False : NonNegativeIntegers
          4 :     0 :  15.0 :  None : False : False : NonNegativeIntegers

  Objectives:
    value : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True : 19.26358596079164

  Constraints:
    const1 : Size=1
        Key  : Lower : Body  : Upper
        None :  None : 100.0 : 100.0
    const2 : Size=4
        Key : Lower : Body : Upper
          1 :  15.0 : 15.0 :  None
          2 :  20.0 : 20.0 :  None
          3 :  50.0 : 50.0 :  None
          4 :   5.0 : 15.0 :  None

4つのタスクに対して、15:20:50:15と割り振ることでパフォーマンスが最大化され、パフォーマンス和は19.26です。さっさと終わる仕事の4番目ですが、プロトタイプが出来てからの伸びしろが良いのでこうなるわけです。

おわり

以上Pythonでの数値最適化でした。皆様も是非ガシガシ最適化して、ついでにブログもアップして、勉強させて下さい。

参考

*1:解説が薄いのは、このライブラリ(記事)を見つけたのが記事を書き終わってからというだけであり、他意はありません

*2:APIが纏められていない点が使いづらいです。出版されている本に載っているようですが、高い…

*3:昨年末に書かれた記事があったので、心配はしていませんでした

Pythonによる超初歩的な金融資産解析(ついでのビットコイン)

概要

資産を増やす金融商品として、投資信託や株、債権なんかがメジャーです。初歩的なポートフォリオ理論では

  • 株などの資産がどのように変動するかは予測することが出来ない
  • 一方で経済は成長するので、全体を長期的に見たらプラスに成長する
  • なので、分散して投資することで安定的な利益を得ることが出来る

といったことが前提とされます。即ち、価格の変動が正規分布に従ってランダムウォークするとみなして、観測された変動の平均と分散を元に、どれだけリスクを取ることが出来るかを踏まえて資産の内訳を作ります。そのため様々な株式や債権を纏めた商品である投資信託は、理論的には優等生な商品と見なすことができます*1

正直な所、金融商品については、記載情報を見ても何を買うべきなのかよく分かりません。一方で、金融では時系列データが結構揃っていて弄ってみるだけでも面白そうです。今回はデータのハンドリングを通して、初歩的な金融資産のバランシングについて解析しました。

注意

  • マサカリ歓迎です。特にこの記事の筆者は金融業界について全くの素人です。橘玲氏の本を読んだ程度の知識しかありません。
  • この記事によって生じたあらゆる事象について、筆者は責任を負いません。
  • この記事は個人の見解と趣味の産物であり、投資を勧めるものでは全くありません。ほんとに。

材料

投資信託基準価格

以下のものをそれぞれのWebサイトからダウンロードしました。

つまり、国内系ファンドはモーニングスター、海外ETFはY! Financeです。国内ファンドはSBIで購入することが出来るものしかありませんから、そちらからもダウンロードできます。しかし自分の環境*2では、csvファイルなのに違う拡張子としてダウンロードされたり、ファイルフォーマットが分析に使いづらいものだったので利用をやめました。モーニングスターはWebサイトが重いという欠点はありますが、ヘッダーがすっきりしていて使いやすいです。

ファンドの選定基準としては、以下のとおりです。

為替価格

海外ETFについては、ドル価格を円価格に修正する必要があります。今回はマネースクウェア・ジャパンからダウンロードしました。より長期のデータソースを知っている人が居たら教えてくれると助かります…。

実際には為替手数料がこの価格に加わる気がしますが、明るくないので割愛しました。

手法

いつものごとくPythonでやります。先ずダウンロードしてきたcsvファイルを読み込んでデータフレームを作ります。Y!から持ってきたデータはいいのですが、モーニングスターのそれはheaderが日本語です。なので、事前にエクセルなんかで英語に書き換えておかないとpandasが怒ります。注意して下さい。

基本的に終値、あるいは調整後終値で計算します。

import pandas as pd
import numpy as np

# USD/JPY比
exchange = pd.read_csv("USDJPY.csv", index_col=0, parse_dates=True)
exchange = exchange[["close"]]

# 日本系投資信託
smbc_astock = pd.read_csv("smbc_astock.csv", index_col=0, parse_dates=True)
smbc_astock.columns = ["smbc_astock"]

hifumi = pd.read_csv("hifumi.csv", index_col=0, parse_dates=True)
hifumi.columns = ["hifumi"]

topix = pd.read_csv("topix.csv", index_col=0, parse_dates=True)
topix.columns = ["topix"]

nikkei = pd.read_csv("nikkei.csv", index_col=0, parse_dates=True)
nikkei.columns = ["nikkei"]

smbc_jcredit = pd.read_csv("smbc_jcredit.csv", index_col=0, parse_dates=True)
smbc_jcredit.columns = ["smbc_jcredit"]

# 海外ETF
# USDで売られているので、為替価格の調整をする
vti = pd.read_csv("VTI.csv", index_col=0, parse_dates=True)
vti = vti[["Adj Close"]]
vti.columns = ["vti"]
vti = pd.DataFrame(vti["vti"] * exchange["close"], columns = ["vti"])
vti = vti.dropna()

vt = pd.read_csv("VT.csv", index_col=0, parse_dates=True)
vt = vt[["Adj Close"]]
vt.columns = ["vt"]
vt = pd.DataFrame(vt["vt"] * exchange["close"], columns=["vt"])
vt = vt.dropna()

vwo = pd.read_csv("VWO.csv", index_col=0, parse_dates=True)
vwo = vwo[["Adj Close"]]
vwo.columns = ["vwo"]
vwo = pd.DataFrame(vwo["vwo"] * exchange["close"], columns=["vwo"])
vwo = vwo.dropna()

# 全部結合する
df = pd.DataFrame([])
df = df.join([
                smbc_astock,
                hifumi,
                topix,
                nikkei,
                smbc_jcredit,
                vti,
                vt,
                vwo,
    ], how="outer")

# 日付情報を初日分から生成       
date = pd.DataFrame(index=pd.date_range(start=df.index[0], end="2018-01-30", period=1))
df = df.join(date, how="outer")
# dataframe用のラベル
order = [
            "smbc_astock",
            "hifumi",
            "topix",
            "nikkei",
            "smbc_jcredit",
            "vti",
            "vt",
            "vwo"
]

df = df[order]

続いて収益率を計算していきます。収益率は購入時の価格と比較した際の、得られる利益のパーセント値です。たとえば100円で購入した資産が110円になったら収益率は10%ですね。

今回はX[年|日]後の収益率を計算し、更に過去Y[年|日]間の収益率から、収益率の統計値を算出します。XとYについては、短期の動きと長期の動きのどちらを見たいのかで、短期: (X=Y=90日)、長期: (X=365年, Y=365年)としました。

今回の計算は最初のデータが有る日から、全日程を計算対象としています。ただし土日祝祭日は証券取引所が開いていないため、NaNになります。とはいえpandasが統計値を計算する時には、これを自動で取り除いてくれます。各商品毎に収益率を計算しているのは、日本とアメリカでNaNになる部分が異なるからです。

統計値については平均、分散、そしてシャープレシオを計算します。シャープレシオは平均/分散で定義される統計値です。高いほど良い金融資産というようにみなされるんだとか。

import matplotlib.pyplot as plt

# 収益率の計算
diff = 90
for i, o in enumerate(order):
    ds = df[o].dropna()
    rds = (ds.iloc[diff:].reset_index(drop=True) - ds.iloc[:-diff].reset_index(drop=True)) / ds.iloc[:-diff].reset_index(drop=True) * 100
    rds.index = ds.index[diff:]
    rdf = pd.DataFrame(rds)
    if i == 0:
        r = rdf
    else:
        r = r.join(rdf, how="outer")
r = r.join(date)


# 統計値の計算
ave_return_list = []
std_return_list = []
sharpe_ratio_list = []
term = 90
for i in range(len(r) - term):
    ave_return_list.append(r.iloc[i:term+i].mean())
    std_return_list.append(r.iloc[i:term+i].std())
    sharpe_ratio_list.append((r.iloc[i:term+i].mean() / r.iloc[i:term+i].std()).to_dict())
ave_return_df = pd.DataFrame(ave_return_list, index=r.index[term:])[order]
std_return_df = pd.DataFrame(std_return_list, index=r.index[term:])[order]
sharpe_ratio_df = pd.DataFrame(sharpe_ratio_list, index=r.index[term:])[order]

# 可視化
fig = plt.figure(figsize=(15,10))

ax = fig.add_subplot(3,2,1)
ave_return_df.plot(title="Long term",
                   xticks=[],
                   ax=ax,
                   legend=False)
ax.set_ylabel("Return")

ax = fig.add_subplot(3,2,2)
ave_return_df.plot(title="Short term",
                   xticks=[],
                   xlim=("2016-10-01", "2017-1-30"),
                   ylim=(-5,15),
                   ax=ax,
                   legend=False)

ax = fig.add_subplot(3,2,3)
std_return_df.plot(xticks=[],
                   ylim=(0,40),
                   ax=ax,
                   legend=False)
ax.set_ylabel("Risk")

ax = fig.add_subplot(3,2,4)
std_return_df.plot(xticks=[],
                   xlim=("2016-10-01", "2017-1-30"),
                   ylim=(0,10),
                   ax=ax,
                   legend=False)

ax = fig.add_subplot(3,2,5)
sharpe_ratio_df.plot(ax=ax,
                     ylim=(-5,20),
                     legend=False)
ax.set_ylabel("Sharpe ratio")

ax = fig.add_subplot(3,2,6)
sharpe_ratio_df.plot(xlim=("2016-10-01", "2017-1-30"),
                     ylim=(0, 4),
                     ax=ax,
                     legend=True)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

plt.show()

最後に、計算した収益率からリスクリターン分布のプロットをします。この図は、横軸に収益率の分散、縦軸に収益率の平均をとったものです。詳細はWikipediaを見てもらうとして、今回挙げた金融商品を一定の割合で保持していた場合に、過去365日間での収益はいかほどだったのかを調べます。

金融資産の割合は、ディリクレ分布からパラメータを適当に変えて生成してみます。最終的に、シャープレシオの順にソートし、1年前に買っていたら今頃儲かっていたという、金融資産の割合を調べます。乱数から統計値を計算する際に、高速化のために逆ソートしてから計算しています。ソートせずに後ろからテーブルを出力しようとすると遅いです。

可視化は最終的に出力したテーブルから出しても良いのですが、そちらだと色が一色になって重なりが見えづらいです。

from scipy.stats import dirichlet
import matplotlib.pyplot as plt
%matplotlib inline

# 収益率の計算
diff = 365
for i, o in enumerate(order):
    ds = df[o].dropna()
    rds = (ds.iloc[diff:].reset_index(drop=True) - ds.iloc[:-diff].reset_index(drop=True)) / ds.iloc[:-diff].reset_index(drop=True) * 100
    rds.index = ds.index[diff:]
    rdf = pd.DataFrame(rds)
    if i == 0:
        r = rdf
    else:
        r = r.join(rdf, how="outer")
r = r.join(date)
[f:id:ajhjhaf:20180203000841p:plain]

# 資産割合の乱数を生成し、その度に可視化
fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(1,1,1)
result = []
term = 365
for alpha in [0.5, 5, 50, 500]:
    for i in range(2000):
        w = dirichlet.rvs(alpha=[alpha / df.shape[1]]*df.shape[1])[0]
        ave_r = r.sort_index(ascending=False).iloc[:term].mean()
        var_r = r.sort_index(ascending=False).iloc[:term].cov()
        e_ret = np.dot(w, ave_r)
        e_ris = np.sqrt(np.dot(np.dot(w, var_r), w.T))
        d = {"weight": w, "return": e_ret, "risk": e_ris, "ratio": e_ret / e_ris}
        result.append(d)
        ax.plot(e_ris, e_ret, "o", ms=3)
ax.set_xlabel("Risk(Standard deviation)")
ax.set_ylabel("Return(Average)")
ax.plot([0, 40], [0, 40])

# 集計
result_df = pd.DataFrame(result)
result_df = result_df.sort_values("ratio", ascending=False)
weight_df = (pd.DataFrame(result_df["weight"].values.tolist(), columns=order) * 100).round()
weight_df.index = result_df.index
result_weight_df = weight_df.join(result_df.sort_values("ratio", ascending=False)[["ratio", "return", "risk"]] / diff * 365)

結果

あまり大した考察は出来ませんが、コメントです。

90日での値動き

f:id:ajhjhaf:20180203063146p:plain

左側がx軸の時系列を長めに取った図、右側が最近の時系列に注目した図です。一番上の段が収益率平均、中段が収益率分散、下段がシャープレシオです。

  • 収益率について、2008年ぐらいにリーマン・ショックがあって株価がめちゃめちゃに落ち込んでいます。しかしそこからは基本的に正のリターンを持っているので、儲かっていますね。ここ最近の収益率も右肩上がりです。

  • こういう図をプロットするまで信じていなかったのですが、たしかに債権はリーマン・ショックのときにも落ち込んでいませんし、資産担保が出来ています。まぁ、一方で全然プラスになっていませんが、安定的だということがわかります。2段目の分散も低いですね。

  • VTとSMBC 全海外は日本が含まれていないこと以外はポートフォリオが似ていますが、実際パフォーマンスも似ています。一方で、ここ最近は日本の株価が上昇傾向だったから、VTの方がパフォーマンスが良かったんですね。

  • VTIはVTよりパフォーマンスが良いと言われることもありますが、ここ最近は分散がかなり高いです。リターンも実際高いんですが、リスクの兼ね合いが大事です。

  • VWOの傾向から、ここ最近は新興国株が安定して上昇傾向だったことがわかります。リターンは高くとも、分散が低く、シャープレシオが凄いことになっていますね。

365日での値動き

f:id:ajhjhaf:20180203063209p:plain

  • ひふみプラスが凄いことになっています。そりゃー人気が出るわけです。ただし、90日のプロットの方を見ればわかりますが、ここ最近のパフォーマンスはやや低めです。

  • 90日の時には殆ど見られませんでしたが、365日だとファンドが始まった時にシャープレシオが安定しないように見えます。基本的に規模が小さすぎて、様子見状態で分散が小さいことが原因と推測します。

  • 短期の場合と異なり、VWOの分散が最も高いです。ただしこちらのほうが、多く言われている状態に一致するようです。その分最近の新興国株式の状態が良かったと推測できます。

  • 儲かっている時には資産がそちらに移動するようで、新興国株と日本債権のシャープレシオが見事に逆相関を描いています。

リスクリターン分布

f:id:ajhjhaf:20180203063240p:plain

sharpe ratio順にソートしたものを上位10個だけ出力します。見やすくするために、1%以下の数値については丸めています。

index smbc_astock hifumi topix nikkei smbc_jcredit vti vt vwo ratio return risk
1386 0 16 0 0 84 0 0 0 5.20125 5.88951 1.13233
313 0 12 0 3 84 0 0 0 5.09683 5.3542 1.0505
1108 0 12 1 0 83 0 4 0 4.97517 5.50391 1.10628
1103 0 6 0 7 86 0 0 0 4.95763 4.2625 0.859785
1302 0 11 0 0 87 1 1 0 4.95401 4.94603 0.998389
1495 0 4 0 8 86 0 2 0 4.93894 3.98458 0.806767
1470 0 12 0 0 88 0 0 0 4.50332 4.88915 1.08568
416 7 0 5 0 87 0 1 0 4.37415 3.30528 0.755639
1355 0 0 0 13 86 0 0 0 4.23473 3.36452 0.794506
  • 効率的フロンティア曲線っぽいものが見えます。が、どこか歪な形ですね。今回は全部で8000個の乱数を生成しているんですが、出し方が完全にランダムなので効率があまり良くありません。シャープレシオが高くなるようにサンプリングを工夫すれば、もっとはっきりと見えるようになると思います。

  • 表を単純に出力しただけでは、リスクをなんとか下げようとしてしまい、債権ばっかになってしまいます。これはこれで理想的なポートフォリオなんですが、リターンは大きくありません。試しにリターンが20%以上のバランスだけ見てみます。

index smbc_astock hifumi topix nikkei smbc_jcredit vti vt vwo ratio return risk
1280 0 73 0 0 26 1 0 0 1.66061 20.6777 12.4519
1945 0 76 0 1 23 0 0 0 1.63506 21.3648 13.0667
217 0 75 0 0 24 0 0 1 1.63401 21.1787 12.9612
169 0 80 0 0 17 0 2 0 1.59045 22.6507 14.2417
176 0 87 0 0 12 1 0 0 1.58291 24.1981 15.2871
1440 0 84 0 0 14 0 0 2 1.5753 23.5152 14.9275
1246 0 87 0 1 12 1 0 0 1.57514 24.2402 15.3892
617 0 89 0 0 10 1 0 0 1.57479 24.7418 15.7112
1264 0 89 0 0 10 0 1 0 1.57079 24.71 15.7309

ひふみは面白いですね。インデックスファンドに勝利してるアクティブファンドって点で凄いと思いました。ただし、これが過去1年間の後ろ向き分析ということには注意して下さい。これから購入しても、同じパフォーマンスが得られるとは限りません。

  • このプロットには、自分のポートフォリオを載せることも出来ます。以下を乱数のプロットが終わった後に記述すればOKです。リバランスするときなんかに、役立てられると思います。
# 保有資産の割合を入力
w = np.array([0.5, 0, 0, 0, 0, 0.35, 0.15, 0])
expected_return = np.dot(w, ave_r)
expected_risk = np.sqrt(np.dot(np.dot(w, var_r), w.T))
sharpe_ratio = expected_return / expected_risk
print("{0}日の期待リターンは{1}%でした".format(diff, expected_return))
print("{0}日の期待リスクは{1}%でした".format(diff, expected_risk))
print("{0}日のシャープレシオは{1}%でした".format(diff, sharpe_ratio))
ax.plot(expected_risk, expected_return, "^", ms=20)

f:id:ajhjhaf:20180203063421p:plain

おまけ

ビットコインについての価格変動も、このスキームで調べてみましょう。今流行っているアレです。

ビットコイン価格

Yahoo! Financeからダウンロードしました。この価格はBTC-USDであることに注意してください。国内の投資信託と合わせるために、為替価格で補正する必要があります。

結果

ビットコインは値動きが激しいことが予測されるので、X=Y=30日での統計値を見ます。

f:id:ajhjhaf:20180203063834p:plain

面白いことに、なんだかんだ言ってここ最近購入していても儲かる可能性はあったようです。ただし分散が30日間で平均してもずいぶん大きいです。なので購入するタイミング次第で物凄い高値を掴みうるリスクがあります。僕が友達と会って「最近ビットコインがアツいけど、もう手遅れだよね」と話したのが11月中頃でしたが、そういえば随分伸びたものです。しかしリスクが凄まじい。場合によっては資産が60%以上丸ごと減る可能性もあるわけです。

続いて直近の30日の履歴から求めたリスクリターン分布のプロットです。

f:id:ajhjhaf:20180203064100p:plain

随分面白い形になりました。超ハイリターンハイリスク資産であるビットコインが入るとこんな感じになるんですね。この場合の最適ポートフォリオは次のとおりです。

index smbc_astock hifumi topix nikkei smbc_jcredit vti vt vwo bc ratio return risk
3856 2 43 0 1 29 0 2 20 2 87.9018 73.4374 10.1646
864 0 0 3 16 2 0 0 75 5 85.8413 94.9392 13.4562
2906 36 23 2 0 17 2 6 12 2 84.7379 68.2586 9.80057
3413 4 21 0 17 0 12 1 42 3 84.1403 90.7351 13.1203
2903 45 17 1 3 7 4 2 19 2 83.465 73.1678 10.6656
2375 52 21 0 3 5 2 0 14 2 83.228 75.7125 11.068
2958 4 52 8 0 8 0 10 15 2 82.6065 90.1786 13.2819
3859 6 15 0 1 1 15 14 45 3 82.3415 89.5655 13.2341
2399 7 0 2 1 6 0 0 79 5 81.6456 97.8794 14.5858

先程までのポートフォリオとは随分異なります。興味深いのはVWOの割合が高くなった点です。ビットコイン新興国株も投機的なところがありますので、この2つを持つことでリスクヘッジすることが出来るんでしょうか笑。仮想通貨界隈では、本格的に市場に出回る前のものをバラバラに買って、†神に祈る†といった投資方法があるようですが、このポートフォリオでもごく少量だけ保持することが示唆された点でも、それは結構合理的なのかもしれません。

とはいえ、仮想通貨での収入は雑収入扱いになるので、最大で50%の所得税がかかります。その点を考慮すると、また違った結果になりそうで面白いですね*3

履歴

*1:ただし、野村の解説にあるように、正規分布に従うランダムウォークは既に否定されています。そもそも変動が正規分布に従うなら、左右対称になっていつまでも経済が成長しません。なので個人的にも歪み正規分布と見なすほうが自然に思えます。「初歩的」と上述したのはそういうことです。

*2:Mac + Safari

*3:計算のどこに取り入れればいいのか分からない