Edwards in his review of Sylla

and with the aid of a calculator, 
I ploughed through Bernoulli's arithmetic 
only to disagree with his answer. He finds 
the expectation to be 4 349/3596 but I find 
4 153/17980 (4.0971 and 4.0085). Bernoulli 
remarks that since the cost of each throw 
is set at 4 "the player's lot is greater 
than that of the peddler' but according to my 
calculation, only by one part in about 500. 
I should be glad to hear from any reader 
who disagrees with my result.


# re: Bernoulli 'error' Significance june 10 2013

nummi = c(120,100,30,24,18,10,6,6,6,5,3,3,3,2,
          1,
          3,3,3,3,4,4,6,8,12,16,24,25,32,180)
          
casus = c(1,16,52,128,245,416,664,976,1369,1776,2204,2560,2893,3088)
Casus = c( casus, 3184, casus[14:1] )

N = sum(Casus) ;  N 
32*31*30*29/(4*3*2*1)                                                                 

e = sum(Casus*nummi)/sum(Casus) ; e

(e -4 )*3596

4 + 349/3596 ; e

cbind(nummi,Casus, nummi*Casus)

4*35960 + 3490 = 147330

E.Num  = 4*17980 + 153
E.Den  = 17980

4*35960 + 306 = 144146
147330 - 144146


sum(Casus*nummi)/sum(Casus)  # just as reported by Bernoulli

####

#  See Bernoulli's Derivation in his addendum to his page 172;
# (available from jh, email below)
#  he only considers the half range, as distribution is symmetric

# This R enumeration, done quite differently, agrees with his..

n=rep(0,32)  # first 3 will remain 0, since Puneta range is 4:32

for(i1 in 1:8) {
	f1 = 4 ;      # 4 possibilites of an i1
	for(i2 in 1:8) {
		f2 = f1 * ( 3 + (i1 != i2) ) # 3/4 possibilites of an i2, depending...
		for (i3 in 1:8) {
		    f3 = f2 * ( 2 + (i3 != i1) + (i3 != i2) )  # etc
		    for (i4 in 1:8) {
		        f4 = f3 * ( 1 + (i4 != i1) + (i4 != i2) + (i4 != i3) ) # etc
			    n[i1+i2+i3+i4] = n[i1+i2+i3+i4] + f4
			}  
		}
	}
}

freq = n[4:32]/24 
cbind(Casus, freq, Casus - freq) ; 

e = sum(nummi*freq)/sum(freq) ; e

(e -4 )*3596

4 + 349/3596 ; e

# The answer Bernoulli gives in the book is 4 + 349/3596, which is 4.0971.

### simulation also agrees with Bernoulli

SIMS=1000000
Total =0 ; Total.sq=0
pocket.values = rep(1:8,4)
INDICES = 1:length(pocket.values)
REPLACE = FALSE
for(sim in 1:SIMS) {
	indices=sample(INDICES,4,replace=REPLACE)
	total.value = sum(pocket.values[indices])
	payout = nummi[total.value-3]
	Total = Total + payout
    Total.sq = Total.sq + payout^2
}
mean= Total/SIMS  # 
var = Total.sq/SIMS - mean^2
mean + c(-1.96,0,1.96)*sqrt(var/SIMS)

# (obtain E of approx 4.5 if allow replacement ie 2 of more in same pocket

########## So, where did Bernoulli and I go wrong ?

#### james.hanley@mcgill.ca