Easy to type

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

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

概要

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

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

 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分布は恒等関数や無限和が入ってきて扱いづらい