주사위 조합의 표준 편차

Question

120까지 합한 주사위 조합 (30)에서 추출 된 숫자의 시퀀스에 대해 stdev를 찾으려고합니다. 파이썬에 매우 익숙하므로이 코드는 콘솔을 고정시킵니다 왜냐하면 숫자는 끝이 없기 때문에 더 작고 효율적인 함수에 모두 맞출 방법을 모르겠습니다. 내가 한 일은 다음과 같습니다 :

• 가능한 모든 조합은 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) 

Answer 1

당신은 6 ³⁰ = 2의 필터링되지 않은 길이의 결과를 요청 * 불가능 ²³ (10)는, 같은 처리합니다.

결합 될 수있는 두 가지 가능성이 있습니다

가

1. 가 포함 더에 생각하고 문제를 사전에 치료, 예를 들어, 합계가 120 인 경우
2. 대신 개의 모든 조합을 샘플링하지 않고 100035의 임의의 쌍만 샘플링하여 대표적인 샘플을 얻으면 충분히 정확하게 표준을 결정할 수 있습니다.

지금, 나는 단지 무력 코드 제공, (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 dices 제품의 표준 편차를 벗어난 경우, 그게 당신의 코드 않습니다. 그렇다면 표준 (1 자리)에 대략 재현 할 수있는 값을 얻기 위해 1 000 000 개의 샘플이 필요합니다. PC를 약 20 초 정도 걸리며 여전히 백만 년 미만입니다 :-D.

같은 것이 있습니까 3.22 * 10 무엇을 찾고 계십니까? 코멘트 후

편집 - 제약 (합 = 120, 총 개수 = 30)에 대입함으로써도 4 실제로 음 숫자의 주파수를 샘플링하는 대신 단지 6 독립 변수를 제공한다. 이 약 1 초에 계산

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.28737023733e + 17로 얻을 수 있습니다 : 내 현재 코드는 다음과 같습니다. 이전의 접근 방식이나이 방법은 버그가 있거나 둘 다 있습니다.

미안 - 그다지 쉽지 않음 : 표본 추출의 가능성은 같지 않습니다. 즉, 여기에서의 문제입니다. 각 샘플은 가능한 조합의 수를 달리하여 무게를 제공하므로 표준 편차를 취하기 전에 고려해야합니다. 위의 코드에서 초안을 작성했습니다.

Answer 2

참고 : 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