やみとものプログラミング日記 やみとものプログラミング日記
TOP 任意の確率密度f(x)に従う連続型確率変数Xのシミュレート方法
任意の確率密度f(x)に従う連続型確率変数Xのシミュレート方法

任意の確率密度f(x)に従う連続型確率変数Xのシミュレート方法

プログラミング 統計
作成日時: 2020年7月11日
更新日時: 2020年7月15日
こんにちは、やみともです。
機械学習を理解したい、Kaggleで強くなりたい、ということで統計を勉強しています。

数学は少し苦手なので、得意なプログラミングと絡めて勉強するのが良いと分かりました。
統計では確率変数を扱いますが、ふと、連続型の確率変数ってプログラムでどうシミュレートするのか気になりました。
今回はその方法を紹介します。

連続型確率変数のシミュレート

言い忘れましたが、言語はPythonで実装します。

今回は例として、$-10 \leq x \leq 10$で定義された$f(x) = \frac{3}{2000}x^2$という確率密度に従う確率変数Xをシミュレートしてみます。
ちなみに$\frac{3}{2000}$という数は$-10 \leq x \leq 10$の範囲で全確率が1になるように計算して出しました。

プログラム

from scipy import stats
import numpy as np

class Quadratic(stats.rv_continuous):
    def _pdf(self, x):
        return (3.0/2000.0) * (x**2)
    
q = Quadratic(a=-10.0, b=10.0)
arr = np.array(q.rvs(size=10))
print(np.round(arr, 1))

科学計算を行うPythonライブラリscipyの統計モジュールstatsを使います。
基本的な流れは、まずstats.rv_continuousという連続確率変数クラスを継承した自前のクラス(上の例ではQuadratic)を用意し、_pdf()メソッドを実装します。
その後、Quadraticクラスをインスタンス化しrvs()メソッドで確率変数の実現値を得るという流れです。

pdfやrvsが何を意味しているのかは次のサイトがまとまっていて分かりやすいです。基本的なscipy.statsモジュールのインターフェースも学べるのでおすすめです。

scipy.stats - scipyの統計関数群のAPI

q = Quadratic(a=-10.0, b=10.0)
の部分について、aとbには確率変数の値が取りうる範囲を設定します。
aが最小値(デフォルトは-∞)、bが最大値(デフォルトが+∞)です。

参考にしたサイト