9

次の Python コードと出力があります。

>>> import numpy as np
>>> s = [12.40265325, -1.3362417499999921, 6.8768662500000062, 25.673127166666703, 19.733372250000002, 21.649556250000003, 7.1676752500000021, -0.85349583333329804, 23.130314250000012, 20.074925250000007, -0.29701574999999281, 17.078694250000012, 3.3652611666666985, 19.491246250000003, -0.76856974999999039, -1.8838917499999965, -6.8547018333333085, 4.5195052500000088, 5.9882702500000136, -9.5889237499999922, 13.98170916666669, -2.3662137499999929, 12.111165250000013, -6.8334957499999902, -21.379336749999993, 8.4651301666666967, 2.5094612500000082, -0.21386274999998989, 5.1226162500000072, 14.283680166666699, -4.3340977499999909, -2.7831607499999933, 8.2339832500000085, -12.841856749999991, -6.4984398333333075, -6.2773697499999912, -13.638411749999996, -15.90088974999999, -8.2505068333333043, -19.616496749999996, -4.4346607499999919, -10.056376749999991, -13.581729833333299, -8.2284047499999957, -4.5957137499999945, -5.3758427499999968, -12.254779833333302, 11.207287250000007, -12.848971749999997, -14.449801749999992, -17.247984749999993, -17.475253833333305]
>>> np.mean(s)
1.3664283380001927e-14
>>> np.std(s)
12.137473069268983
>>> (s - np.mean(s)) / np.std(s)
array([ 1.02184806, -0.11009225,  0.56658138,  2.1151954 , ...

これを R で実行すると、結果が一致しません。

> options(digits=16)
> s = c(12.40265325, -1.3362417499999921, 6.8768662500000062, 25.673127166666703, 19.733372250000002, 21.649556250000003, 7.1676752500000021, -0.85349583333329804, 23.130314250000012, 20.074925250000007, -0.29701574999999281, 17.078694250000012, 3.3652611666666985, 19.491246250000003, -0.76856974999999039, -1.8838917499999965, -6.8547018333333085, 4.5195052500000088, 5.9882702500000136, -9.5889237499999922, 13.98170916666669, -2.3662137499999929, 12.111165250000013, -6.8334957499999902, -21.379336749999993, 8.4651301666666967, 2.5094612500000082, -0.21386274999998989, 5.1226162500000072, 14.283680166666699, -4.3340977499999909, -2.7831607499999933, 8.2339832500000085, -12.841856749999991, -6.4984398333333075, -6.2773697499999912, -13.638411749999996, -15.90088974999999, -8.2505068333333043, -19.616496749999996, -4.4346607499999919, -10.056376749999991, -13.581729833333299, -8.2284047499999957, -4.5957137499999945, -5.3758427499999968, -12.254779833333302, 11.207287250000007, -12.848971749999997, -14.449801749999992, -17.247984749999993, -17.475253833333305)
> mean(s)
[1] 1.243449787580175e-14
> sd(s)
[1] 12.25589024484334
> (s - mean(s)) / sd(s)
 [1]  1.01197489551755737 -0.10902853430514588  2.09475824715945480  0.56110703609584245 ...

違いはかなり小さいことはわかっていますが、これは私のアプリケーションでは少し問題です。また、R の結果は Stata の結果とも一致することに注意してください。

注: Python 2.7.2、NumpPy 1.6.1、R 2.15.2 GUI 1.53 Leopard ビルド 64 ビット (6335) を使用しています。

4

3 に答える 3

9

sこれは、元の投稿で与えられたリストを使用して、プレーンPythonを使用してその一部に光を当てます:

>>> import math
>>> sum(s) / len(s)
1.3664283380001927e-14
>>> math.fsum(s) / len(s)
1.2434497875801753e-14

最初の出力は を再現np.mean()し、2 番目の出力は R を再現します(R コードが使用されていれば、それらは同一になるとmean()確信しています)。options(digits=17)

Python での違いは、sum()"左から右へ" 追加することで、加算のたびに丸め誤差が生じますが、math.fsum()概念的には無限精度の合計を計算し、最後に 1 つの丸めの総計を使用して、無限精度の合計を最も近い表現可能なものに置き換えます。倍精度数。

Dollars to donuts は、R もそれを行っていると述べています。それは、@John が、数値の順序に関係なく R が同じ平均を返すと報告する理由を説明しますs(無限精度の合計は、加数の順序にまったく影響されません)。

それで終わりだとは思いませんが。Rはおそらくstd devの計算にもより優れた数値的方法を使用しています.数値誤差が小さいという意味では「より良い」が、計算により多くの時間がかかるという意味ではおそらく「悪い」.

PEP 450 - 「統計モジュールを標準ライブラリに追加する」が最近 Python に受け入れられたことに注意してください。これにより、これらのものの高品質な (数値的に) 実装が標準ライブラリに追加されます。もちろん、numpyこれらも使用するかどうかはユーザー次第です。

標準開発について

平均はどのように計算されても 0 に近く、 の数値sはまったく 0 に近くないため、計算された平均の差はほとんど無関係です。これを証明するために、無限精度の計算を行う構成要素を次に示します (これも単純な Python です)。

from fractions import Fraction
def sumsq(xs):
    fs = [Fraction(x) for x in xs]
    mean = sum(fs) / len(fs)
    return sum((f - mean)**2 for f in fs)

これを使用して、母集団とサンプルの標準偏差の非常に高品質の (そして非常に遅い!) 推定値を生成できます。

>>> ss = sumsq(s)
>>> ss  # exact result:  no rounding errors so far!
Fraction(606931231449932225838747590566767, 79228162514264337593543950336)
>>> from math import sqrt
>>> sqrt(ss / len(s))  # population sdev with 2 roundings
12.137473069268983 
>>> sqrt(ss / (len(s) - 1))     # sample sdev with 2 roundings
12.255890244843338

だから - 驚き、驚き ;-) -母集団の標準偏差に対する可能な限り最良np.std(s)の二重近似を計算し、Rはサンプル標準偏差に対する可能な限り最良の二重近似を計算しました。sd()

したがって、この特定のケースでは、計算された平均値間の数値の差は厄介者であり、元の数値に比べて平均値が小さかったため、標準偏差を計算するほぼすべての方法で良好な数値結果が得られます。

ここでの本当の違いは、R がデフォルトで分母 (サンプル sdev) を使用するのに対し、デフォルトで分母 (母集団 sdev) をnp使用することだけです。nn-1

于 2013-10-05T00:26:58.233 に答える
5

64 ビットの精度は 2e-16 程度しかないことに注意してください。これらの数値を合計すると、平均と同様に合計が 0 に非常に近いことがわかります。したがって、問題はその精度に関係している可能性があります。参照する各関数は、最初に数値を合計する必要があります。ということで、もとに戻りました。

RReduce('+', s)では、python 関数と同じ合計が得られますsum。R と Python 内では、実際にはまったく同じように合計されます。ただし、R の関数meanとsum関数は、より正確な方法を使用して計算を行います。numpyで行われているのと同じ方法でR内ですべての数学を行うと、同じになります。

使用している Python 計算について懸念する理由があります。あなたが使用している R コードは、実際には物事をより適切に処理しています。試す：

# R
sum(s)
sum(s * 10000) / 10000
Reduce('+', s)
Reduce('+', s*10000)/10000

# python (numpy is the same here)
sum(s)
sum(s * 10000) / 10000

両方のsum合計が同じであるため、R ではスケーリングが適切に処理されます。ただし、R も python も逐次合計法を使用してそれを処理することはできません。あなたが試みることができるもう一つのことは、数字をスクランブルすることです. コードは提供しませんがsum、R では一貫して同じ値を与えますが、RReduceとsumpython では注文に応じて異なる値を与えます

それで、あなたは何をしますか？精度が非常に高いことを受け入れ、0に近い値を0として扱う必要があることをお勧めします.偏差。合計から得られる平均誤差は、分散を開始してから爆発するだけです。おそらく、そのような数値が同じでなければならない正確な理由についての詳細情報は、より正確なアドバイスを得るのに役立ちます。

同一であることがすべて重要な場合に機能する代替手段があります。R の組み込み関数を使用しないでください。それらは高品質すぎて、でこぼこした統計の問題を浮き彫りにしています。合計を示したように平均値と標準偏差をロールするとReduce、結果は同じになります。しかし、あなたがしようとしているのは、R を遅くし、精度を下げることです。このオプションをまったく回避できる場合は、そうしてください。例えば：

npMean <- function(x) Reduce('+', x)/length(x)
npMean(s)
npSD <- function(x) {m <- npMean(x); sqrt( Reduce('+', (x - m)^2)/(length(x)) )}
npSD(s)

正確にpython平均と（間違った）numpy SDを提供します。それらは機能しますが、R の根性を回避するのが難しい場合があり、物事が正確になりすぎます。もちろん、numpy 関数を置き換えて Python コードをより正確にする Python 関数を見つけることができれば、それはさらに良いことです。

于 2013-10-04T23:02:58.320 に答える