やりたいこと

P10、最頻値（モード）、P90が分かっているとき、その条件にフィットする三角分布を生成する

三角分布は、最小値a、最頻値b、最大値cの三つのパラメータで表現される連続分布の一種です。¹
モデルが単純で直感的に理解しやすいため、Oil&Gas業界のリスク評価などでよく用いられています。²

ここでは、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()

出力

0.10000000062440002 0.8999999953027843

参考にしたサイト

最小値a、最大値b、最頻値cと定義されている場合もあります ↩
たとえば「石油上流業界のポートフォリオ戦略」など ↩