やりたいこと
- P10、最頻値(モード)、P90が分かっているとき、その条件にフィットする三角分布を生成する
三角分布は、最小値a、最頻値b、最大値cの三つのパラメータで表現される連続分布の一種です。1
モデルが単純で直感的に理解しやすいため、Oil&Gas業界のリスク評価などでよく用いられています。2
ここでは、P10、最頻値、P90が既知の標本データがあるとき、これにフィットする三角分布を計算する方法を紹介します。
SciPyのcurve_fitを用いて最小二乗法によるフィッティングを行います。
コード
次のような三角分布のクラスtriangular_fitを定義します。
triangular_dist.py
import numpy as np
from scipy.stats import triang
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
#三角分布のクラスを定義
class triangular_fit:
def __init__(self,p10,mode,p90):
#与えられたP10,最頻値、P90にフィットする三角分布を求める
self.b = mode
array_x = [p10,p90]
array_y = [0.1,0.9]
#curve_fitを用いて、累積分布関数がP10,P90にフィットする三角分布のパラメータa,cを計算
param, cov = curve_fit(self.cdf_fit, array_x, array_y,bounds=[(-np.inf,p90),(p10,np.inf)])
self.a = param[0]
self.c = param[1]
def cdf_fit(self,x,a,c):
#累積分布関数を出力(フィッティング計算用)
#scipyの三角分布の関数用にパラメータを変換
c_tri=(self.b-a)/(c-a)
loc_tri=a
scale_tri=c-a
return triang.cdf(x,c_tri,loc=loc_tri,scale=scale_tri)
def cdf(self,x):
#累積分布関数を出力
#scipyの三角分布の関数用にパラメータを変換
c_tri=(self.b-self.a)/(self.c-self.a)
loc_tri=self.a
scale_tri=self.c-self.a
return triang.cdf(x,c_tri,loc=loc_tri,scale=scale_tri)
def pdf(self,x):
#確率密度関数を出力
#scipyの三角分布の関数用にパラメータを変換
c_tri=(self.b-self.a)/(self.c-self.a)
loc_tri=self.a
scale_tri=self.c-self.a
return triang.pdf(x,c_tri,loc=loc_tri,scale=scale_tri)
def ppf(self,q):
#パーセントポイント関数(CDFの逆関数)を出力
#scipyの三角分布の関数用にパラメータを変換
c_tri=(self.b-self.a)/(self.c-self.a)
loc_tri=self.a
scale_tri=self.c-self.a
return triang.ppf(q,c_tri,loc=loc_tri,scale=scale_tri)
ここでは、P10=0、Mode=5、P90=15を満たす三角分布を計算してみます。
p10 = 0
mode = 5
p90 = 15
#与えられたP10,mode,P90にフィットする三角形分布
test = triangular_fit(p10,mode,p90)
#チェック用:P10、P90の累積分布値がそれぞれ0.1,0.9に近い値になっているか確認。
print(test.cdf(p10),test.cdf(p90))
#チェック用作図:確率密度関数を出力
x = np.linspace(0-p90/2,p90*2,1000)
plt.plot(x,test.pdf(x))
plt.show()
出力
参考にしたサイト
-
最小値a、最大値b、最頻値cと定義されている場合もあります ↩
-
たとえば「石油上流業界のポートフォリオ戦略」など ↩