ダイスの出目によるポイントの合計値の出現率を求めたい

Question

前提

1. 10面ダイスをX個振る
2. それぞれのダイスについて、出目が1だったら2、10だったら-1、それ以外の場合はY以下なら1、Yより大きいなら0をポイントとする
3. Yが1～9のそれぞれの場合において、ポイントの合計の出現率を求めたい

例えば、

X=1, Y=5の場合

ポイント	出現率
-1	10%
0	40%
1	40%
2	10%

X=3, Y=7の場合

ポイント	出現率
-3	0.1%
-2	0.6%
-1	3%
0	8.3%
1	19.2%
2	26.4%
3	29.1%
4	11.4%
5	1.8%
6	0.1%

となります。
Xを1～30程度、Yを1～9の範囲での全ての結果を求めたいと思っています。

試したこと

まずは高速に動作するであろうC++で愚直に実装してみました。

#include <iostream>
#include <array>
#include <cmath>

int main()
{
	// ダイスの出目によるポイントのテーブル
	const char point_table[10][10] = {
		{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
		{2, 0, 0, 0, 0, 0, 0, 0, 0, -1},
		{2, 1, 0, 0, 0, 0, 0, 0, 0, -1},
		{2, 1, 1, 0, 0, 0, 0, 0, 0, -1},
		{2, 1, 1, 1, 0, 0, 0, 0, 0, -1},
		{2, 1, 1, 1, 1, 0, 0, 0, 0, -1},
		{2, 1, 1, 1, 1, 1, 0, 0, 0, -1},
		{2, 1, 1, 1, 1, 1, 1, 0, 0, -1},
		{2, 1, 1, 1, 1, 1, 1, 1, 0, -1},
		{2, 1, 1, 1, 1, 1, 1, 1, 1, -1},
	};

	const int X = 8;

	// X個のダイスを振った時の全ての組み合わせにおいて、合計ポイントが出現する回数を記録する配列
	// ポイントの合計は-X～2Xの範囲で出現するのでそれが収まる範囲
	// point_sum[X]が合計ポイント0の出現回数で、X-1なら合計ポイント-1、X+2なら合計ポイント2の出現回数を指す
	std::array<int, X*3+1> point_sum = {0};

	for (int Y = 1; Y <= 9; Y++) {
		for (size_t i = 0; i < point_sum.size(); i++) {
			point_sum[i] = 0;
		}
		for (int i = 0; i < 10; i++) {
			int p1 = point_table[Y][i];
			for (int j = 0; j < 10; j++) {
				int p2 = point_table[Y][j];
				for (int k = 0; k < 10; k++) {
					int p3 = point_table[Y][k];
					for (int l = 0; l < 10; l++) {
						int p4 = point_table[Y][l];
						for (int m = 0; m < 10; m++) {
							int p5 = point_table[Y][m];
							for (int n = 0; n < 10; n++) {
								int p6 = point_table[Y][n];
								for (int o = 0; o < 10; o++) {
									int p7 = point_table[Y][o];
									for (int p = 0; p < 10; p++) {
										int p8 = point_table[Y][p];
										// p1～p8はそのダイスによるポイント
										// その合計は-X～2Xの範囲で出現するので、配列に収まるようそれにXを加算したものを配列のインデックスにしている
										point_sum[p1+p2+p3+p4+p5+p6+p7+p8+X]++;
									}
								}
							}
						}
					}
				}
			}
		}
		for (size_t i = 0; i < point_sum.size(); i++) {
			std::cout << "X = " << X << " / Y = " << Y << " / " << (int)i-X << " = " << point_sum[i]/std::pow(10, X-2) << "%" << std::endl;
		}
	}

	return 0;
}

発生している問題

このプログラムはX=8の場合の、出現する全てのダイス目をカウントする事で確率を求めています。
この愚直な方法でも確かに確率は求められるのですが、当然ですがダイスの数が1つ増えるごとに実行時間がおよそ10倍になり、こちらの環境ではX=10の段階で計算終了まで20分ほどかかってしまいました。

また、あるXの時に計算した内容を保持しておき、X+1の計算に流用する事で計算時間を減らせないかとも検討しました。
上記のX=8の例であれば、i, j, k, l, m, n, o, pの全ての組み合わせにおいてp1+p2+p3+p4+p5+p6+p7+p8の値を記録しておき、X=9の計算の際にその値を利用する方法も試したのですが、僅かに実行時間は減らせるものの結局はXが増えるたびに指数的に実行時間が伸びるのであまり意味はなく、メモリ的な限界もあり難しいと判断しました。

実現したいこと

ダイスの数が30個程度（できれば、少なくとも20個）までの結果は求めたいと思っているのですが、この方法で求められるのはせいぜい12個程度までで、それ以上は実行時間が非現実的です。
なにか良い方法はありますでしょうか？
目的は確率を求める事なので、実行時間は現実的な範囲であれば長くなっても問題ありませんし、プログラミング以外の手段でそれが求められるのならそれでも構いません。
言語、ライブラリ等も問いません。

Answer 1

とりあえず頭空っぽにして考えると，
１個のダイスで加算されるポイントのパターンは４パターンしかないので，
X個のダイスを考えたときには，4^X パターンの結果を総当たりすれば済む，すなわち

ダイスの数が1つ増えるごとに実行時間がおよそ10倍になり

の「10倍」のところを，4倍にはできるのではないでしょうか．

要は，「出目の確率が偏った４面ダイス」を扱うのだと考えて良かろう，と．
（この「４面ダイス」の各目の出る確率はYで定まる）

---

で，愚直にX回の繰り返しをするとしたら，
各回の時点で同じ得点のパターンがあればそいつらを統合してパターンを減らしてやればいいのでは．

例えば，X=3のとき，愚直にやると

1. １個目のダイスの結果が４パターン
2. ２個目のダイスでさらに４パターンだから，ここまでの組み合わせは 4*4=16 パターン
3. ３個目のダイスでさらに４パターンだから，組み合わせは 16*4 = 64 パターン

となるかもしれないけども，実際には2.の時点で取り得る総ポイントのパターンは -2 ～ 4 の7パターンしかないので，ここでデータをまとめてしまえば，3.での組み合わせを 7*4 = 28 パターンに減らせる．

---

愚直に実装してみた．
（「早い」かどうかは謎だけど20分とかはかからないのではないかと）

//Yによって定まる，１個のサイコロから得られるポイントの確率データを作る作業用関数．
//（めんどくさいのでYはまともな値が入ってくる前提）
void SetupPointRate( int Y, double (&Rates)[4] )
{
	Rates[0] = 0.1;  //-1 point の確率
	Rates[2] = (Y-1) * 0.1;	//1 point
	Rates[3] = 0.1;	// 2 point

	Rates[1] = 0.8 - Rates[2];	//0 point
}

//処理実装
struct Data
{
public:
	Data() : m_MinPoint(0), m_P{1.0} {}
public:
	//サイコロ１回分のデータ更新作業
	void Update( const double (&Rates)[4] )
	{
		const std::vector<double> PrevP = m_P;
		const int PrevMinPoint = m_MinPoint;
		--m_MinPoint;
		m_P.assign( m_P.size()+3, 0 );
		
		for( int iPrev=0; iPrev<PrevP.size(); ++iPrev )
		{
			int Point = PrevMinPoint + iPrev;
			for( int iRate=0; iRate<4; ++iRate )
			{
				int deltaPoint = iRate - 1;
				m_P[ Point+deltaPoint-m_MinPoint ] += ( PrevP[iPrev] * Rates[iRate] );
			}
		}
	}

	//結果表示
	void Show() const
	{
		int Point = m_MinPoint;
		for( const auto &P : m_P ){	std::cout << Point++ << " : " << P*100 << "%\n";	}
	}
private:
	int m_MinPoint;	//現時点での最低ポイント
	std::vector<double> m_P;	//各ポイントの確率
};

//main
int main(void)
{
	const int X=3;
	const int Y=7;

	double Rates[4];
	SetupPointRate( Y, Rates );
	Data D;
	for( int i=0; i<X; ++i ){	D.Update( Rates );	}
	D.Show();

	return 0;
}

Comments

[fana]

ダイス１個増やす度に総ポイントのパターンが３パターン増えるという点は
質問コードで既に用いられていますね．

[fana]

とりあえず回答に書いたコードで，質問文に例示されている２例
・X=1, Y=5
・X=3, Y=7
については，結果が質問文記載の表と一致することを確認した．

実行時間に関しては実行環境依存となろうが，
X=30, Y=7 とかでも「一瞬」で終わるので，Yのパターン分繰り返しても大した時間は要しないのではないかと思う．
（結果があってるかどうかの検算まではしていないが．面倒すぎるので^^）

[angeltime]

回答ありがとうございました。
考え方だけでなくコードまでつけて頂き、非常に分かりやすかったです。
自分のコードで求めたX=10までの結果とも等しかったです。
とても助かりました！

Answer 2

サイコロ１回の確率分布はYによって最初に決まるので、あとはX回累積していけばよいでしょう。
そのさい、同じポイントならまとめます。
以下、Pythonで愚直に実装しても一瞬です。

# サイコロ一回の確率分布
def get_p(y):
    d = {} # キー=ポイント, 値=確率
    d[2] = 0.1
    d[-1] = 0.1
    d[1] = 0.1*(y-1)
    d[0] = 0.1*(9-y)
    return d

def calc(X,Y):
    cs = {0:1} # 初期はポイント０で確率は１
    ps = get_p(Y)
    for x in range(X):
        cn = {} # 累積確率
        for kc,vc in cs.items():
            # 次のサイコロの結果
            for kp,vp in ps.items():
                k = kc + kp
                v = vc * vp
                if k not in cn:
                    cn[k] = 0
                cn[k] += v # 同じポイントなら重ね合わせ
        cs = cn

    print(f'X={X}, Y={Y}')
    for k in sorted(cs.keys()):
        print(f'{k:2d} = {cs[k]*100}%')
    total = sum(cs.values())
    print(f'total={total*100}%')

calc(1,5)
"""
X=1, Y=5
-1 = 10.0%
 0 = 40.0%
 1 = 40.0%
 2 = 10.0%
total=100.0%
"""

calc(3,7)
"""
X=3, Y=7
-3 = 0.10000000000000002%
-2 = 0.6000000000000002%
-1 = 3.0000000000000004%
 0 = 8.300000000000002%
 1 = 19.200000000000006%
 2 = 26.400000000000006%
 3 = 29.10000000000001%
 4 = 11.400000000000004%
 5 = 1.8000000000000005%
 6 = 0.10000000000000002%
total=100.00000000000003%
"""

calc(30,7)
"""
X=30, Y=7
-30 = 1.0000000000000017e-28%
-29 = 6.00000000000001e-27%
-28 = 1.920000000000003e-25%
-27 = 4.295000000000007e-24%
-26 = 7.482000000000012e-23%
-25 = 1.0752852000000018e-21%
-24 = 1.3212907500000025e-20%
-23 = 1.4221402800000024e-19%
-22 = 1.3642783200000026e-18%
-21 = 1.1818021246000021e-17%
-20 = 9.338154924960017e-17%
-19 = 6.784952085000013e-16%
-18 = 4.562878534900509e-15%
-17 = 2.855474942272205e-14%
-16 = 1.6704082999252832e-13%
-15 = 9.169100281793844e-13%
-14 = 4.73808035180431e-12%
-13 = 2.3113356502258042e-11%
-12 = 1.0669826300269449e-10%
-11 = 4.670853808253032e-10%
-10 = 1.942555783975026e-09%
-9 = 7.687447347558051e-09%
-8 = 2.898865920221126e-08%
-7 = 1.0428957093227026e-07%
-6 = 3.583292565890607e-07%
-5 = 1.1769385860713239e-06%
-4 = 3.6983040975543256e-06%
-3 = 1.1125676547788542e-05%
-2 = 3.2061039402944635e-05%
-1 = 8.854550263163403e-05%
 0 = 0.00023445760643797724%
 1 = 0.0005953950945956934%
 2 = 0.0014504017071929436%
 3 = 0.003389836242944168%
 4 = 0.0076016699979969266%
 5 = 0.01635609743971916%
 6 = 0.03376445030939296%
 7 = 0.0668631008684591%
 8 = 0.12698902083101876%
 9 = 0.23124522828867378%
10 = 0.40359702466316716%
11 = 0.6748387159590861%
12 = 1.0804451869333935%
13 = 1.6553669852579709%
14 = 2.4253553557748746%
15 = 3.395530033498209%
16 = 4.5384882925340495%
17 = 5.785796255731776%
18 = 7.02739157052915%
19 = 8.122438699582188%
20 = 8.92217522980884%
21 = 9.30087041063324%
22 = 9.186806600477004%
23 = 8.583292848029442%
24 = 7.57165170052293%
25 = 6.293775796550219%
26 = 4.91915190657864%
27 = 3.6069620346182942%
28 = 2.475273254616503%
29 = 1.5857826981064225%
30 = 0.9459604951271544%
31 = 0.5240296824383176%
32 = 0.26886166661489375%
33 = 0.1274204777116115%
34 = 0.055637019724470695%
35 = 0.022326780386830936%
36 = 0.008214996935576495%
37 = 0.0027653630129821235%
38 = 0.0008498934857256249%
39 = 0.0002380112317889737%
40 = 6.062265577455562e-05%
41 = 1.4017432498517185e-05%
42 = 2.9367646696764194e-06%
43 = 5.563607991254611e-07%
44 = 9.509638554506276e-08%
45 = 1.4628530592320676e-08%
46 = 2.0193348162987405e-09%
47 = 2.4929366639816986e-10%
48 = 2.7412972803885247e-11%
49 = 2.6720008404564667e-12%
50 = 2.2950647489833977e-13%
51 = 1.7246053638910043e-14%
52 = 1.1235949518480028e-15%
53 = 6.274837938000016e-17%
54 = 2.9597725235000067e-18%
55 = 1.1562150360000026e-19%
56 = 3.6400800000000086e-21%
57 = 8.874300000000019e-23%
58 = 1.5720000000000033e-24%
59 = 1.8000000000000032e-26%
60 = 1.0000000000000017e-28%
total=100.00000000000028%
"""

Comments

[angeltime]

回答ありがとうございます。
今回は先に回答頂いた方をベストアンサーにしましたが、コードの内容としてはそちらとほぼ同じような考え方かと思いますので、こちらも非常に参考になる回答だと思います。
ありがとうございました。

Answer 3

配列を数値じゃなくて確率にして、一個前の値ごとの確率を掛ければ結構短縮できそう、
例えばX=4,X=7のとき、2~-1となる確率はpoint_table={0.1, 0.6, 0.2, 0.1}で
例のX=3,Y=7の結果からｰ3~6となる確率がX3={0.001,0.006,0.03,0.083,0.192,0.264,0.291,0.114,0.018,0.001}だと判明しているのなら
X3の1番目からpoint_tableを順にかけた値をX4の配列に加えるのを最後まで繰り返せば早そう。
X4={0.0001,0.0006,0.0002,0.0001,0,0,0,0,0,0,0,0,0}
X4={0.0001,0.0006+0.0006,0.0002+0.0036,0.0001+0.0012,0+0.006,0,0,0,0,0,0,0,0}
-3のときの4つ目のダイスの結果とｰ2の時の4つ目のダイスの結果は別事象なので足し合わせられるはず。

Comments

[angeltime]

回答ありがとうございます。
今回は先に回答頂いた方をベストアンサーにしましたが、考え方のヒントを教えてくださり勉強になります。
ありがとうございました。

Answer 4

openmpを使ってみてはどうでしょうか

Comments

[angeltime]

回答ありがとうございます。
今回は解決できるコードを提示してくださった方が居たため、そちらの方法で行いたいと思います。
openmpについても今後のために調べておきたいと思います。