1

計算したい

ここに画像の説明を入力

nまでの値を1000000可能な限り正確に。ここにいくつかのサンプルコードがあります。

from __future__ import division
from scipy.misc import comb
def M(n):
    return sum(comb(n,k,exact=True)*(1/n)*(1-k/n)**(2*n-k)*(k/n)**(k-1) for k in xrange(1,n+1))
for i in xrange(1,1000000,100):
    print i,M(i)

最初の問題はOverflowError: long int too large to convert、n = 1101. これは、comb(n,k,exact=True)が大きすぎて float に変換できないためです。ただし、最終結果は常に 前後の数値になり0.159ます。

大きな中間値を持つ合計を計算する方法で関連する質問をしましたが、この質問は主に 3 つの理由で異なります。

  1. 計算したい数式が異なるため、さまざまな問題が発生します。
  2. 私が示した例に見られるように、exact=True を使用するために以前に提案された解決策は、ここでは役に立ちません。浮動小数点除算を実行する必要があるため、comm の独自の実装をコーディングしてもうまくいきません。
  3. 以前よりもはるかに大きな値の答えを計算する必要があるため、新しい問題が発生します。何らかの巧妙な方法で合計をコーディングしないと、それはできないと思います。

クラッシュしない解決策は、使用することです

from fractions import Fraction
def M2(n):
    return sum(comb(n,k,exact=True)*Fraction(1,n)*(1-Fraction(k,n))**(2*n-k)*Fraction(k,n)**(k-1) for k in xrange(1,n+1))
for i in xrange(1,1000000,100):
    print i, M2(i)*1.0

残念ながら、現在は非常に遅いためn=1101、妥当な時間内に回答が得られません。

したがって、2 番目の問題は、大きなn.

4

3 に答える 3

4

乗算、除算、累乗をそれぞれ加算、減算、乗算に置き換える対数変換を使用して、各被加数を計算できます。

def summand(n,k):
    lk=log(k)
    ln=log(n)
    a=(lk-ln)*(k-1)
    b=(log(n-k)-ln)*(2*n-k)
    c=-ln
    d=sum(log(x) for x in xrange(n-k+1,n+1))-sum(log(x) for x in xrange(1,k+1))
    return exp(a+b+c+d)

def M(n):
    return sum(summand(n,k) for k in xrange(1,n))

k=n の場合、被加数はゼロになるため、対数が未定義になるため計算しないことに注意してください。

于 2013-04-17T13:12:04.300 に答える
3

You can use gmpy2. It has arbitrary precision floating point arithmetic with large exponent bounds.

from __future__ import division
from gmpy2 import comb,mpfr,fsum

def M(n):
    return fsum(comb(n,k)*(mpfr(1)/n)*(mpfr(1)-mpfr(k)/n)**(mpfr(2)*n-k)*(mpfr(k)/n)**(k-1) for k in xrange(1,n+1))
for i in xrange(1,1000000,100):
    print i,M(i)

Here is an excerpt of the output:

2001 0.15857490038127975
2101 0.15857582611615381
2201 0.15857666768820194
2301 0.15857743607577454
2401 0.15857814042739268
2501 0.15857878842787806
2601 0.15857938657957615

Disclaimer: I maintain gmpy2.

于 2013-04-17T13:18:05.670 に答える
2

かなり残忍な方法は、すべての要因を計算し、結果が約 1.0 にとどまるように掛け算することです (Python 3.x):

def M(n):
    return sum(summand(n, k) for k in range(1, n + 1))

def f1(n, k):
    for i in range(k - 1):
        yield k
    for i in range(k):
        yield n - i

def f2(n, k):
    for i in range(k - 1):
        yield 1 / n
    for i in range(2 * n - k):
        yield 1 - k / n
    yield 1 / n
    for i in range(2, k + 1):
        yield 1 / i

def summand(n, k):
    result = 1.0
    factors1 = f1(n, k)
    factors2 = f2(n, k)
    while True:
        empty1 = False
        for factor in factors1:
            result *= factor
            if result > 1:
                break
        else:
            empty1 = True
        for factor in factors2:
            result *= factor
            if result < 1:
                break
        else:
            if empty1:
                break
    return result

M(1101)が得られますが、数0.15855899364641846秒かかります。M(2000)約 14 秒かかります0.15857489065619598。

（最適化できると思います。）

于 2013-04-17T12:31:54.777 に答える