2015-12-06 33 views
5

Özellikle, bir toplamını hesaplamak R kullanıyorumR neden burada toplamı yanlış yapıyor?

(cevap R hala bu bu konuda bilgisi olan birinden bunun nedenleriyle ilgili ayrıntıları duymak isteriz, küçük ondalık sayıları depolamak için ikili sayıları nasıl kullandığını gerçekten ise) Bu bir: çıktıya bakarsanız

r = 21 #Number of times the loop executes 
n = 52 
sum = 0 #Running total 
for(k in 0:r){ 
    sum = sum + (1-(choose(2^(k), n)*(factorial(n)/(2^(n*k))))) 
    print(sum) 
} 

, fark edeceksiniz:

enter image description here

İşte benim Ar kod

[1] 1 
[1] 2 
... 
[1] 11.71419 
[1] 11.71923 
[1] 11.72176 
[1] 12.72176 
[1] 13.72176 

19. yinelemeden sonra neden 1 artırılıyor?

Bu görev için daha uygun olabilecek başka serbest hesaplamalı motorlar var mı?

+2

Çünkü 1'den sonraki terim, k = 19 ve NaN (underflow) sonra 25'ten sonra 0'dır. Toplamın ne yapmasını bekliyorsunuz? Sonlu bir şeye yakınlaşmak mı? –

+2

Maple'ı deneyebilirsiniz, toplamı yaklaşık 11.72428812 olarak bulur. R açısından, bir olasılık rasgele değişken için Monte Carlo simülasyonu yapmaktır. – Julius

+3

Sadece kayıt için, döngü yazıldığı gibi, döngü r + 1 kez (ek açıklamada belirtildiği gibi r kere değil) yürütülür –

cevap

5

Tüm hesaplamalar burada genellikle fakülte ve yükselme yetkileri için uygun olmayan kayan nokta sayıları kullanmaktadır (çünkü bu değerler hızlı hesaplanır ve bu sayede hesaplama daha az doğru olur).

> factorial(52)/(2^(52*19)) 
[1] 3.083278e-230 
> factorial(52)/(2^(52*20)) 
[1] 0 

ve Pari-GP ile karşılaştırmak:

, deneyin

? \p 256 
realprecision = 269 significant digits (256 digits displayed) 
? precision((52!)/2^(52*19) + . , 24) 
%1 = 3.08327794368108826958435659488643289724 E-230 
? precision((52!)/2^(52*20) + . , 24) 
%2 = 6.8462523287873017654100595727727116496 E-246 
5

Eğer yaşıyorsanız taşma/Yetersizlik sorunu aşmanın bir yolunu arıyorsanız, Eğer (ara hesaplamalar için makul büyüklüklerini tutmak için) günlükleri kullanın ve ardından sonunda exponentiate edebilirsiniz:

options(digits=20) 

for(k in 0:r){ 
    sum = sum + (1 - exp(lchoose(2^k, n) + log(factorial(n)) - k*n*log(2))) 
    print(paste0(k,": ",sum)) 
} 

[1] "0: 1" 
[1] "1: 2" 
[1] "2: 3" 
... 
[1] "19: 11.7217600143979" 
[1] "20: 11.7230238079842" 
[1] "21: 11.7236558993777" 

Bunun doğru olduğunu kontrol etmek için, Mathematica'da orijinal toplamı (günlük çekmeden) çalıştırdım ve aynı sonuçları 12 ondalık basamağa getirdim.

R numaralı sorunu çözebilirseniz de, bir bilgisayar cebir sistemiyle (sembolik matematik ve kesin hesaplamalar yapmanıza olanak verir) gitmek istiyorsanız, Sage ücretsiz ve açık bir kaynaktır.