0

我正在尝试查找从骰子（30）组合中提取的数字序列的 stdev，总和为 120。我对 Python 很陌生，所以这段代码使控制台冻结，因为数字是无穷无尽的，而我不是确定如何将它们全部放入一个更小、更高效的函数中。我所做的是：

  • 找到 30 个骰子的所有可能组合；
  • 过滤后的组合总计为 120；
  • 将结果列表中列表中的所有项目相乘；
  • 尝试提取标准偏差。

这是代码：

import itertools
import numpy

dice = [1,2,3,4,5,6]
subset = itertools.product(dice, repeat = 30)

result = []
for x in subset:
    if sum(x) == 120:
        result.append(x)

my_result = numpy.product(result, axis = 1).tolist()
std = numpy.std(my_result)

print(std)
4

2 回答 2

1

请注意D(X^2) = E(X^2) - E(X)^2，您可以通过以下等式解析地解决此问题。

f[i][N] = sum(k*f[i-1][N-k])        (1<=k<=6)
g[i][N] = sum(k^2*g[i-1][N-k])
h[i][N] = sum(h[i-1][N-k])

f[1][k] = k ( 1<=k<=6)
g[1][k] = k^2 ( 1<=k<=6)
h[1][k] = 1 ( 1<=k<=6)

示例实现：

import numpy as np

Nmax = 120
nmax = 30
min_value = 1
max_value = 6
f = np.zeros((nmax+1, Nmax+1), dtype ='object')
g = np.zeros((nmax+1, Nmax+1), dtype ='object') # the intermediate results will be really huge, to keep them accurate we have to utilize python big-int
h = np.zeros((nmax+1, Nmax+1), dtype ='object')
for i in range(min_value, max_value+1):
    f[1][i] = i
    g[1][i] = i**2
    h[1][i] = 1

for i in range(2, nmax+1):
    for N in range(1, Nmax+1):
        f[i][N] = 0
        g[i][N] = 0
        h[i][N] = 0
        for k in range(min_value, max_value+1):
            f[i][N] += k*f[i-1][N-k]
            g[i][N] += (k**2)*g[i-1][N-k]
            h[i][N] += h[i-1][N-k]

result = np.sqrt(float(g[nmax][Nmax]) / h[nmax][Nmax] - (float(f[nmax][Nmax]) / h[nmax][Nmax]) ** 2)
# result = 32128174994365296.0
于 2016-10-30T00:41:38.053 回答
0

您要求未过滤长度为 6 30 = 2*10 23的结果，因此无法处理。

有两种可能性可以组合：

  1. 包括更多思考来预处理问题，例如关于如何仅对总和为 120 的样本进行抽样。
  2. 而是进行蒙特卡罗模拟，即不要对所有组合进行采样，而仅对 1000 个随机组合进行采样以获得具有代表性的样本，从而确定足够准确的标准。

现在，我只应用（2），给出蛮力代码：

N = 30 # number of dices
M = 100000 # number of samples
S = 120 # required sum

result = [[random.randint(1,6) for _ in xrange(N)] for _ in xrange(M)]
result = [s for s in result if sum(s) == S]

现在，该结果应该与您使用之前的结果相当numpy.product……虽然那部分我无法遵循……

好的，如果您在 30 个骰子乘积的标准差之后出局，这就是您的代码所做的。然后我需要 1 000 000 个样本来获得 std（1 位）的大致可重现值 - 我的 PC 大约需要 20 秒，仍然远少于 100 万年:-D。

像 3.22*10 16这样的数字是你要找的吗？

评论后编辑： 嗯，对数字的频率进行采样只会给出 6 个自变量 - 实际上甚至 4 个，通过替换约束（总和 = 120，总数 = 30）。我当前的代码如下所示：

def p2(b, s):
    return 2**b * 3**s[0] * 4**s[1] * 5**s[2] * 6**s[3]

hits = range(31)
subset = itertools.product(hits, repeat=4) # only 3,4,5,6 frequencies
product = []
permutations = []
for s in subset:
    b = 90 - (2*s[0] + 3*s[1] + 4*s[2] + 5*s[3]) # 2 frequency
    a = 30 - (b + sum(s)) # 1 frequency
    if 0 <= b <= 30 and 0 <= a <= 30:
        product.append(p2(b, s))
        permutations.append(1) # TODO: Replace 1 with possible permutations
print numpy.std(product)  # TODO: calculate std manually, considering permutations

这大约需要 1 秒计算，但令人困惑的是我得到的结果是 1.28737023733e+17。我以前的方法或这个方法有一个错误 - 或两者兼而有之。

抱歉——没那么容易：抽样的概率不同——这就是问题所在。每个样本都有不同数量的可能组合，给出了它的权重，在计算标准偏差之前必须考虑这一点。我已经在上面的代码中起草了这一点。

于 2016-10-29T21:42:59.590 回答