Heres a monte carlo simulator written in R. Personally I like tripling the bet each time, transferring out winnings once they reach a certain level then starting again. The thing with martingales is you are bound to lose eventually if you follow it strictly:
#install.packages("Rlab") # Run First Time
require(Rlab) #used for bernouilli distribution function rbern()
####Settings
#General
live.plot=F #Watch Simulation Live (runs slower)
max.bets<-1000 #Number of Bets to End Simulation
iterations<-1000 #Number of Simulations to Run
fee<-.0005 #Transaction Fee
#Bet Strategy and Satoshi Dice options
Win.Odds<-.5 #Odds of Winning
Price.Multiplier<-1.957 #Multiple of Bet Payed on Win
max.bet.size<-500 #Maximum Allowed Size of Bet
start.wallet<-10 #Starting Funds
lose.multiplier<-2 #Multiple of losing bet size to use for following bet
#End Settings
##Calulate Bet Sizes and Payouts
bet.num<-1:20
bet.size<-.01*lose.multiplier^(bet.num-1)
bet.size[which(bet.size>max.bet.size)]<-max.bet.size
win.size<-Price.Multiplier*.01*lose.multiplier^(bet.num-1)-fee
win.size[which(win.size>Price.Multiplier*max.bet.size)]<-Price.Multiplier*max.bet.size-fee
###Highest Density Interval Calculator Function (for the charts)
get.HDI<-function(sampleVec,credMass){
sortedPts = sort( sampleVec )
ciIdxInc = floor( credMass * length( sortedPts ) )
nCIs = length( sortedPts ) - ciIdxInc
ciWidth = rep( 0 , nCIs )
for ( i in 1:nCIs ) {
ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
}
HDImin = sortedPts[ which.min( ciWidth ) ]
HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
HDIlim = c( HDImin , HDImax, credMass )
return( HDIlim )
}
####
######Run Simulation
if(live.plot==T){
dev.new()
}
out2=NULL
pb<-txtProgressBar(min = 0, max = iterations, initial = 0,style = 3)
for(j in 1:iterations){
wallet<-start.wallet
i<-1
t<-1
run.result=NULL
while(wallet>0 & t<=max.bets){
result<-rbern(1,Win.Odds)
if(result==0){
wallet<-wallet-bet.size[i]
color="Red"
i<-i+1
}
if(result==1){
wallet<-wallet-bet.size[i]+win.size[i]
color="Green"
i<-1
}
if(wallet>0){
t<-t+1
run.result<-rbind(run.result,cbind(j,t,wallet,i))
}
if(live.plot==T){
plot(run.result[,2],run.result[,3],
xlab="Bet Number", ylab="Coins in Wallet",
col=color,
main=paste("Simulation #",j)
)
}
}
out2<-rbind(out2,run.result)
setTxtProgressBar(pb, j)
}
close(pb)
#End Simulation
####Get Results
winning.iterations<-out2[which(out2[,2]==max.bets),1]
perc.wins<-100*length(which(out2[,2]==max.bets))/iterations
win.wallets<-out2[which(out2[,2]==max.bets),3]
hdi<-get.HDI(win.wallets,.95)
win.mode<-density(win.wallets)$x[
which(density(win.wallets)$y==
max(density(win.wallets)$y)
)]
fail.timepoints=matrix(nrow=(iterations-length(win.wallets)),ncol=1)
r<-1
for(i in 1:iterations){
temp<-out2[which(out2[,1]==i),]
if(max(temp[,2]) fail.timepoints[r]<-max(temp[,2])
r<-r+1
}
}
####Plot Results
dev.new()
layout(matrix(c(1,1,2,3),nrow=2,ncol=2,byrow=2))
#Plot Simulation Results
plot(0,0,type="n",
xlab="Bet Number", ylab="Coins in Wallet",
xlim=c(0,max(out2[,2])),ylim=c(0,max(out2[,3])),
main=c(paste("Starting Funds=", start.wallet, " Win Odds=", Win.Odds),
paste("Bet Multiplier After Loss= ",lose.multiplier,"x",sep=""),
paste("Maximum # of Bets=",max.bets, " # of Simulations=",iterations))
)
for(i in 1:iterations){
if(i %in% winning.iterations){
temp<-out2[which(out2[,1]==i),]
lines(temp[,2],temp[,3],col=rgb(0,1,0,.1))
}else{
temp<-out2[which(out2[,1]==i),]
lines(temp[,2],temp[,3],col=rgb(1,0,0,.1))
}
}
#Plot Distribution of non-zero Winnings
hist.counts<-hist(win.wallets, col="Green",
xlab="Final Wallet Amount",
breaks=seq(0,(max(win.wallets)+max(win.wallets)/10),by=max(win.wallets)/10),
main=c(paste("Percent of Winning Simulations=",perc.wins, "%"),
paste("Mode=",round(win.mode,2), " Mean=",round(mean(win.wallets),2))
),
sub=paste(round(100*hdi[3],3), "%", "HDI (Lower,Upper):",round(hdi[1],1),",",round(hdi[2],1))
)$count
rect(hdi[1],0,hdi[2],.025*max(hist.counts),col="Black")
#Plot Distribution of when Bankruptcy Occurred
hist(fail.timepoints, col="Red",
breaks=seq(0,(max(fail.timepoints)+max(fail.timepoints)/10),by=max(fail.timepoints)/10),
xlab="Bet Number at Bankruptcy",
main="Bet Number at Bankruptcy")
