Go to Qiita Advent Calendar Top
LoginSignup
0
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@sakaimaging(Yusuke Sakai)in立教大学

連続ポアソン分布を描く：ガンマ関数での拡張とスプライン・ガウス近似の比較

Posted at

はじめに

ポアソン分布は、観測回数などのカウントデータを扱うときによく登場します。
その定義は次のように表されます。

P(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}

ここで：

  • $\lambda$：単位時間（または領域）あたりの平均発生回数
  • $k$：観測されたイベントの回数（0, 1, 2, … のような離散的な整数値）

を表します。

ただし、この式は $k$ が整数に限定されているため、もし $k$ を連続的に扱いたいケース、たとえば勾配を求めたい場合や最適化問題に応用したい場合などでは、そのままでは使えません。

そこで登場するのが、階乗 $k!$ をガンマ関数で拡張し、$k$ を実数として扱えるようにした「連続ポアソン分布（Continuous Poisson distribution）」です（e.g., Kim et al. 2010, Ilienko 2013）。

P(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{\Gamma(k+1)}, \quad k \ge 0

この形では $k$ を連続変数として扱えるため、ポアソン分布をなめらかに近似した連続関数として利用できます。
実際に、X線天文学の解析でも（e.g., Suzuki et al. 2025）、ポアソン尤度を連続化した手法が応用されています。

本記事では、この連続ポアソン分布を Pythonで描画して、通常の離散ポアソン分布と比べてどの程度滑らかに補間されるのかを可視化します。
さらに、スプライン補間やガウス近似と比較し、連続化の違いがどのように現れるのかを確認していきます。

本記事の実装コードは Google Colab こちら からも確認できます。

階乗を連続化する

まず、ポアソン分布の連続化の出発点となる「ガンマ関数」を見てみましょう。
階乗は整数にしか定義されませんが、ガンマ関数を使うと実数にまで滑らかに拡張できます。

下の図では、青い点が $k!$、オレンジの線が $\Gamma(k+1)$ を表しています。
整数点で両者が完全に一致しており、ガンマ関数が階乗の自然な連続拡張であることがわかります。

image.png

コードはこちら 
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.special import factorial

plt.rcParams["font.size"] = 13

# 比較用データ
x_int = np.arange(0, 8)
x_cont = np.linspace(0, 7.1, 300)

# 計算
fact_int = factorial(x_int)
gamma_cont = gamma(x_cont + 1)

# プロット
plt.figure(figsize=(8, 5))
plt.plot(x_cont, gamma_cont, label="Γ(k+1)", color="C1")
plt.scatter(x_int, fact_int, label="k!", color="C0", zorder=5)
plt.title("factorial vs gamma")
plt.xlabel("k")
plt.ylabel("value")
plt.legend()
plt.grid(alpha=0.3)
plt.show()

この「$k! \to \Gamma(k+1)$」という置き換えこそが、連続ポアソン分布を導く基本的なアイデアです。

離散ポアソンをなめらかに

次に、ポアソン分布を連続的に拡張してみましょう。
$k!$ を $\Gamma(k+1)$ に置き換えることで、$k$ を実数として扱えるようになります。

結果を描画すると、青い点（離散ポアソン）がオレンジの滑らかな曲線（連続ポアソン）で自然に補間されているのが分かります。

image.png

コードはこちら 
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.special import factorial

plt.rcParams["font.size"] = 13

# --- 定義 ---
def continuous_poisson(k, lam):
    return (lam**k * np.exp(-lam)) / gamma(k+1)

def discrete_poisson(k, lam):
    return (lam**k * np.exp(-lam)) / factorial(k)

# --- パラメータ設定 ---
lam = 5
k_cont = np.linspace(0, 15, 300)
k_disc = np.arange(0, 16).astype(int)

# --- 計算 ---
p_cont = continuous_poisson(k_cont, lam)
p_disc = [discrete_poisson(k, lam) for k in k_disc]

# --- プロット ---
plt.figure(figsize=(8, 5))
plt.plot(k_cont, p_cont, label="Continuous Poisson", color="C1", lw=2)
plt.scatter(k_disc, p_disc, color="C0", zorder=5, label="Discrete Poisson")
plt.title(f"Continuous vs Discrete Poisson (λ={lam})")
plt.xlabel("k")
plt.ylabel("P(k; λ)")
plt.legend()
plt.grid(alpha=0.3)
plt.show()

このように、連続ポアソン分布は離散的な確率質量関数を滑らかに補間した形になっており、
勾配計算や最適化アルゴリズムで扱える「連続的なポアソン尤度」として利用できます。

他の連続化手法との比較

最後に、連続ポアソン分布を他の一般的な連続化手法と比較してみます。
ここでは次の2つを取り上げます。

  • スプライン補間（Cubic spline）
    離散点を通るように滑らかに結ぶ古典的な補間手法（詳しくはこちら）。
  • ガウス近似（Normal approximation）
    $\lambda$ を平均・分散とする正規分布で近似する方法（詳しくはこちら）。

image.png

コードはこちら 
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.special import factorial
from scipy.interpolate import interp1d
from scipy.stats import norm

plt.rcParams["font.size"] = 13

# --- 定義 ---
def continuous_poisson(k, lam):
    return (lam**k * np.exp(-lam)) / gamma(k+1)

def discrete_poisson(k, lam):
    return (lam**k * np.exp(-lam)) / factorial(k)

# --- パラメータ ---
lam = 50
sigma = np.sqrt(lam)
k_min = max(0, lam - 3 * sigma - 1)
k_max = lam + 3 * sigma + 1
k_cont = np.linspace(k_min, k_max, 300)
k_disc = np.arange(np.floor(k_min), np.ceil(k_max) + 1)

# --- 計算 ---
p_cont = continuous_poisson(k_cont, lam)
p_disc = discrete_poisson(k_disc, lam)

# --- スプライン補間 & ガウス近似 ---
spline = interp1d(k_disc, p_disc, kind="cubic")
p_spline = spline(k_cont)
p_gauss = norm.pdf(k_cont, loc=lam, scale=np.sqrt(lam))

# --- 差分 ---
diff_spline = p_spline - p_cont
diff_gauss  = p_gauss  - p_cont

# --- プロット ---
fig, axes = plt.subplots(2, 1, figsize=(8, 8), sharex=True)
plt.subplots_adjust(hspace=0.05)

# 上段：分布比較
axes[0].scatter(k_disc, p_disc, color="C0", zorder=5, label="Discrete Poisson")
axes[0].plot(k_cont, p_cont, "C1-", lw=2, label="Continuous Poisson (Γ ext.)")
axes[0].plot(k_cont, p_spline, "C2--", lw=1.5, label="Cubic spline interpolation")
axes[0].plot(k_cont, p_gauss, "C3-.", lw=2, label="Gaussian approximation")
axes[0].set_ylabel("P(k; λ)")
axes[0].set_title(f"Comparison of continuous extensions (λ={lam})")
axes[0].legend(loc="lower center", fontsize="small")
axes[0].grid(alpha=0.3)

# 下段：差分
axes[1].plot(k_cont, diff_spline, "C2--", lw=1.5, label="Spline - Continuous")
axes[1].plot(k_cont, diff_gauss, "C3-.", lw=2, label="Gaussian - Continuous")
axes[1].axhline(0, color="k", lw=1, alpha=0.7)
axes[1].set_xlabel("k")
axes[1].set_ylabel("Difference")
axes[1].set_yscale("symlog", linthresh=1e-7)
axes[1].set_xlim(lam - 3 * sigma, lam + 3 * sigma)
axes[1].set_ylim(-5e-2, 5e-2)
axes[1].legend(loc="upper right", fontsize="small")
axes[1].grid(alpha=0.3)

plt.show()

上段では分布の形状を、下段では連続ポアソンとの差分を示しています。

結果からわかること：

  • スプライン補間は整数点で完全に一致し、全体として連続ポアソンに非常に近い。
    ただし、点と点の間ではわずかにズレが生じる。
  • ガウス近似は $\lambda$ の前後で非対称な歪みが見られる。

λを変化させてみる

さらに、$\lambda$ を 1〜60 の範囲で動かしながら、
$\pm 3\sigma$ の範囲における分布の変化をアニメーションで確認してみます。

continuous_poisson_lambda_variation.gif

コードはこちら 
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.special import factorial
from scipy.interpolate import interp1d
from scipy.stats import norm
from matplotlib.animation import FuncAnimation

plt.rcParams["font.size"] = 13
plt.rcParams["figure.figsize"] = (8, 8)

# --- 定義 ---
def continuous_poisson(k, lam):
    return (lam**k * np.exp(-lam)) / gamma(k+1)

def discrete_poisson(k, lam):
    return (lam**k * np.exp(-lam)) / factorial(k)

# --- λ の範囲 ---
lams = np.arange(1, 61, 1)

# --- Figure 準備 ---
fig, axes = plt.subplots(2, 1, figsize=(8, 8), sharex=True)
plt.subplots_adjust(hspace=0.05)

def update(frame):
    lam = lams[frame]
    sigma = np.sqrt(lam)
    k_min = max(0, lam - 3 * sigma - 1)
    k_max = lam + 3 * sigma + 1
    k_cont = np.linspace(k_min, k_max, 300)
    k_disc = np.arange(np.floor(k_min), np.ceil(k_max) + 1)

    # --- 計算 ---
    p_cont = continuous_poisson(k_cont, lam)
    p_disc = discrete_poisson(k_disc, lam)

    # --- スプライン補間 & ガウス近似 ---
    spline = interp1d(k_disc, p_disc, kind="cubic", fill_value="extrapolate")
    p_spline = spline(k_cont)
    p_gauss = norm.pdf(k_cont, loc=lam, scale=np.sqrt(lam))

    # --- 差分 ---
    diff_spline = p_spline - p_cont
    diff_gauss  = p_gauss  - p_cont

    # --- 軸初期化 ---
    axes[0].cla()
    axes[1].cla()

    # --- 上段：分布比較 ---
    axes[0].scatter(k_disc, p_disc, color="C0", zorder=5, s=12, label="Discrete Poisson")
    axes[0].plot(k_cont, p_cont, "C1-", lw=2, label="Continuous Poisson (Γ ext.)")
    axes[0].plot(k_cont, p_spline, "C2--", lw=1.5, label="Cubic spline interpolation")
    axes[0].plot(k_cont, p_gauss, "C3-.", lw=2, label="Gaussian approximation")
    axes[0].set_ylabel("P(k; λ)")
    axes[0].set_title(f"Comparison of continuous extensions (λ={lam})")
    axes[0].legend(loc="lower center", fontsize="small")
    axes[0].grid(alpha=0.3)
    axes[0].set_xlim(k_min, k_max)
    axes[0].set_ylim(0, max(np.max(p_disc), np.max(p_cont), np.max(p_gauss)) * 1.3)

    # --- 下段：差分 ---
    axes[1].plot(k_cont, diff_spline, "C2--", lw=1.5, label="Spline - Continuous")
    axes[1].plot(k_cont, diff_gauss, "C3-.", lw=2, label="Gaussian - Continuous")
    axes[1].axhline(0, color="k", lw=1, alpha=0.7)
    axes[1].set_xlabel("k")
    axes[1].set_ylabel("Difference")
    axes[1].set_yscale("symlog", linthresh=1e-7)
    axes[1].set_xlim(lam - 3 * sigma, lam + 3 * sigma)
    axes[1].set_ylim(-5e-2, 5e-2)
    axes[1].legend(loc="upper right", fontsize="small")
    axes[1].grid(alpha=0.3)

# --- アニメーション生成 ---
ani = FuncAnimation(fig, update, frames=len(lams), interval=300, blit=False)

# GIFで保存（数十秒かかります）
ani.save("continuous_poisson_lambda_variation.gif", writer="pillow", fps=5, dpi=60)
plt.close()

$\lambda$ が大きくなるにつれて、ポアソン分布は正規分布に近づいていくことが分かります。
また、離散ポアソンをスプライン補間したものも、連続分布へとなめらかに近づいていく様子が興味深いです。

おわりに

ポアソン分布をガンマ関数で連続化すると、
整数しか扱えなかった世界が“なめらか”につながり、
勾配や傾きといった連続的な情報を使った計算にも応用できるようになります。

もちろん、これは少し概念的には変な話です。
ポアソンは本来「離散的」な確率分布であり、
連続化にはそれ相応のリスクもあります。
たとえば、連続ポアソンは確率密度分布ではないため、
そのまま確率として解釈することはできません。

それでも、勾配を使った最適化や確率的シミュレーションのように、
連続的な変化を扱う計算の場面では、
この“連続ポアソン”的な考え方が役に立つことがあります。
その一例として、筆者が最近個人的に関心を持っている
CIM（Coherent Ising Machine） のシミュレータ
（CIMを方程式から眺める記事）も、
離散スピンを連続的な関数で近似して扱うという点で、
発想としてどこか似ていますね。

離散的なものを連続的に“ゆるめて考える”──
結局のところ、どれだけ美しくそのあいだを埋められるかが鍵なのだと思います。
理論が描く近似でも、データが導くスプラインでも、
そのあいだの埋め方にはさまざまな流儀があって、
その上で解かれる最適化問題にもまた多様な方法があり、
まさに最適化するための最適化があるのだと、しみじみ感じます。

0
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?