モンテカルロ法での最大値

Question

やりたいこと

Pythonを利用して、モンテカルロ法の計算を行っています。

考えている問題点

コードは下記に記した通り、各パラメータに分布を与えて数値を発生させたうえで計算する方式を考えてます。

当初はそれらしい結果が出てきたので上手くいっていると思っていたのですが、比較すべき統計値との乖離が大きく、結果を分析してみたところ現状の計算で出てくる最大値(out[6]:577256)の結果と各パラメータの最大値で無理やり計算した結果(Out[7]:625780)の乖離が原因では無いかと考えています。
現状のコードではリストに収納された数値同士の組み合わせでの計算なので、必ずしもすべてのパラメータが最大値のケースは発生しないことが原因なのではと思っています。

考えている対応策

対応策としては、1.無理やり全通り計算させるような方式に変更するか、2.乱数の生成数（現状10000）と実際の計算トライアルの数を異なるものとすれば(例：乱数の生成を1000通り、本計算の回数を10万通り)、改善されるのではと思っていますが、Python初心者でしてうまいコードの書き方が判りません。

アドバイスよろしくお願い致します。

###Python
コード
#Monte Carlo
In [1]:import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('darkgrid')

#Properties
In [2]:trials = 10000

A = np.random.uniform(0.1, 0.3, trials)
B = np.random.uniform(500, 700, trials)
C = np.random.triangular(10, 20, 30, trials)
D = 100
In [3]:x = A * B * C * D

x = x.astype(int)
count = np.bincount(x)
mode = np.argmax(count)
In [4]:fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(x, bins=18, density=False, rwidth=0.6, Color="steelblue")
ax.set_title('Test')
ax.set_xlabel('x')
ax.set_ylabel('Frequency')

Out[4]:Text(0, 0.5, 'Frequency')
In [5]:print('P90: {:.0f} '.format(sorted(x)[int(trials * 0.1)]))
print('P50: {:.0f} '.format(sorted(x)[int(trials * 0.5)]))
print('P10: {:.0f} '.format(sorted(x)[int(trials * 0.9)]))
print('Mean: {:.0f} '.format(int(np.mean(x))))
print('Stdv: {:.0f} '.format(int(np.std(x))))
print('Mode: {:.0f} '.format(int(np.argmax(count))))

P90: 132027
P50: 229595
P10: 366067
Mean: 241177
Stdv: 90505
Mode: 199554

In [6]:max(x)
Out[6]:577256

In [7]:x2 = max(A) * max(B) * max(C) * D
x2
Out[7]:625780.9249227662

Comments

[ForestSeo]

コードは
```Python
コード
```
このような形で投稿してください

[Nyann]

ありがとうございます。ご指摘頂いた点修正させていただきました。
よろしくお願いします。

[m.ts10806]

分からない方はこちら参考に。
https://teratail.com/questions/238564

[Nyann]

ありがとうございます。次回投稿する際にはこちらを確認したうえで投稿させていただきます。

[m.ts10806]

いえ、次回では確実に忘れるので、今。

Answer 1

コードの問題というよりも確率的起こりうるやむを得ないことかと思います。

統計量の計算対象であるxは、A×BとCによって生成されていると分解できます。（Dは定数なので、分布に影響しないので省略します）
このうちCは、使用している分布とパラメータから正規分布と似た特性（平均と最頻値が同じ・平均を中心に分布が左右対称）を持つのでモンテカルロシミュレーションの結果も平均に収束します。一方、A×Bは台形のような分布と取ります。質問に記載の1万件だとわかりにくいのですが、10万件で試行すると80～150が一様分布を示します。つまり、A×Bについては80～150の範囲については安定性に欠けるため、計算上の平均である120にうまく収束しないことが予測されます。そのため、A×B×Cもそれなりの誤差が生じてしまうのではないでしょうか。

対策としては、試行回数を増やすことが一番いいのではないでしょうか。つまり、A×Bの分布の一様分布が真に一様分布と同じようになれば、不安定性がなくなるので2400にうまく収束すると思います。

Comments

[Nyann]

ありがとうございます。試行回数を増やしても1000回以上では大きな差は得られませんでした。Pythonコードに大きな問題が無いことが判ったので、比較対象にしているシミュレーション値の算出に問題がある可能性が高いと考えています。先方に確認中ですので判明しましたら情報更新させていただきます。

Answer 2

現状の結果のどの部分に問題を感じているのか不明なのですが、現状の分布でシミュレートすると、99.9パーセンタイル値が 539000 あたりに、99.95パーセンタイル値でさえも 553000 あたりとなるので、たかだか 10000個のサンプルで最大値が 577256 となるのは妥当なのではないでしょうか。

---

【コメントを受けてソースコードの変更】

# Monte Carlo
import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns 
import statistics
sns.set_style('darkgrid')

# Properties
trials = 10000000
A = np.random.uniform(0.1, 0.3, trials)
B = np.random.uniform(500, 700, trials)
C = np.random.triangular(10, 20, 30, trials) # 三角分布
D = 100

x = A * B * C * D

x = x.astype(int)

fig, ax = plt.subplots()
ax.hist(x, bins=18, density=False, rwidth=0.6,  Color="steelblue")
ax.set_title('Test')
ax.set_xlabel('x')
ax.set_ylabel('Frequency')
# plt.show()

mean = np.mean(x)
std = np.std(x)
mode, _ = stats.mode(x) 
p0,p10,p50,p90,p100 = np.percentile(x, q=[0,10,50,90,100])

print(f'P90: {p90:.0f}')
print(f'P50: {p50:.0f}')
print(f'P10: {p10:.0f}')
print(f'Mean: {mean:.0f}')
print(f'Stdv: {std:.0f}')
print(f'Mode: {mode[0]:.0f}')

# 基本統計量だけならば下記のコードでも良い（MODEはないけど）
print(pd.Series(x).describe())
# scipy.stats を使って基本統計量を出した場合は歪度と尖度も出る
print(stats.describe(x))

Comments

[Nyann]

ありがとうございます。確かに統計的にはそうなると思います。
投稿に使用している数字・計算は模擬なのですが、実際のデータではMean、Median共に基準に対して9%程度の乖離が出ているので投稿致しました。
それほどの開きがあると手法(計算式等)自体に問題がある可能性も高いと思うのでそちらも検証するつもりではおりますが、まずは現状のPythonコードに問題無いか、また投稿させて頂いた対応策ではどのような結果が得られるかも確認したいと考えております。

[magichan]

ソースコードに関しては大きな問題は無いかと思いますが、モードやパーセンタイルの計算などにおいて時間がかかる処理を行っております（サンプル数を大きく取ると顕著になるかと思います。）
パーセンタイルの計算はnumpy.percentile() を、モードの計算は scipy.stats.mode() を使うことをおすすめします。
また基本統計量は pandas.Series.describe() や scipy.stats.describe() などでも導けますので、こちらを使ってもよいかとおもいます。
＊回答にサンプルを追記しておきました

[magichan]

あと実際のデータとシミュレート上のデータの平均の乖離の件ですが、こちらは９％という値はあまり意味を持たないのかと思います。シミュレートにて論理値が得られるのであれば、統計的にその母集団の平均値との実際の標本平均との間で平均値の検定を行い、その差が有意かどうかで判断するべきではないでしょうか。

[Nyann]

コードのアドバイスありがとうございました。確かにこちらの方がスマートなので活用させていただきます。
今回比較対象としている結果もシミュレーション値(おそらくCrystal Ball)でして、基本的には全く同じ計算を行っていることとなるため、その割には乖離が大きいと感じています。コードに問題が無いとすれば、適用している数式(若しくは入力値)に問題があるか、先方のミスの可能性があると思いますのでその点先方へ確認しているところです。
今回のようなケース、もしも各計算式に入力されている数値を出力することが出来れば検証しやすいと思うのですが、その場合どのようなコードを入力すればよいかアドバイス戴けますと大変助かります。よろしくお願いいたします。

[magichan]

コードで１点気になる点があるとしたら
x = x.astype(int)
の部分でしょうか。この処理を入れている理由がわからないので何とも言えないのですが、計算で得られるデータは浮動小数なのに対し全てのデータを整数に丸めていますので、この処理が平均値などに影響を与えることは間違いありません。
> 各計算式に入力されている数値を出力することが出来れば検証しやすい
具体的に何の値をどのように表示したいのでしょうか？

[Nyann]

x = x.astype(int)ですが、削除しました。以前のコードでModeを計算するために入れたコードでして、修正頂いたコードでは不要と判りました。
出力については変数の出力を意味しておりました。下記のコードの通り、変数のdfを作成しdf.to_csvで出力を試みているのですが、出力されません。コードの読み込みはうまくいきます。他のネットで出ているサンプルコードですと問題なく出力されます。
データ数の問題かと思い、trialを100に減らしてみましたが改善されませんでした。

--------
df = pd.DataFrame([A, B, C],index= ["A","B","C"])
df = df.T
df.to_csv("XXXX.csv")

[Nyann]

失礼しました。よく見たところJupyterのHomeからはcsvファイルが確認できるので出力はされているようです。確認不足で申し訳ありません。

[magichan]

x = A * B * C * D
の次の行にでも
df = pd.DataFrame({'A':A, 'B':B, 'C':C, 'D':D, 'x':x})
df.to_csv('XXXX.csv')
のようにしてみてください

[Nyann]

ありがとうございます、見事に配列されました。transformを使わなくてもこれで出来てしまうんですね。