In [1]:
using Gadfly
In [2]:
function down(n,r)
    return prod([n-r+1:n])
end
Out[2]:
down (generic function with 1 method)
In [3]:
n = 365
vals = Float64[down(big(n),r)/big(n)^r for r=1:100]
plot(x=[1:100],y=1-vals,Geom.line)
Out[3]:
x -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 -1 0 1 2 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 y
In [4]:
function Stirling(n)
    return n^n*exp(-n)*sqrt(2π*n) 
end
Out[4]:
Stirling (generic function with 1 method)
In [5]:
# Comparing values of factorials to Stirling's approximations 
hcat([Stirling(n) for n=1:10],[factorial(n) for n=1:10])
Out[5]:
10x2 Array{Float64,2}:
      0.922137       1.0     
      1.919          2.0     
      5.83621        6.0     
     23.5062        24.0     
    118.019        120.0     
    710.078        720.0     
   4980.4         5040.0     
  39902.4        40320.0     
 359537.0       362880.0     
      3.5987e6       3.6288e6
In [6]:
function randompermutation(n::Integer)
    numbers = [1:n]
    result = Int64[]
    for i=1:n
        push!(result,numbers[rand(1:length(numbers))])
        numbers = filter( x-> x != result[end], numbers)
    end
    return result
end
Out[6]:
randompermutation (generic function with 1 method)
In [7]:
randompermutation(10)
Out[7]:
10-element Array{Int64,1}:
  9
  8
  1
  6
  5
  2
  4
  7
 10
  3
In [8]:
# count the number of fixed points, i.e., numbers k which are 
# in position k
mean([sum(randompermutation(10) .== [1:10]) == 0 for _=1:10_000])
Out[8]:
0.3782
In [9]:
1/e
Out[9]:
0.36787944117144233
In [10]:
# counts the number of records in a permutation, i.e., numbers which 
# are bigger than all the preceding ones 
function records(permutation)
    r = 0
    result = 0
    for i=1:length(permutation)
        if permutation[i] > r
            result += 1
            r = permutation[i]
        end
    end
    return result
end
Out[10]:
records (generic function with 1 method)
In [11]:
mean([records(randompermutation(200)) for _=1:1000])
Out[11]:
5.868
In [12]:
function BinomialProbabilities(n,p)
    return Float64[binomial(big(n),r)*p^r*(1-p)^(n-r) for r=0:n] 
end

function PowerCurve(n,p,m)
    # Returns the probability of getting m or more successes in n Bernoulli
    # trials of probability p
    return sum(BinomialProbabilities(n,p)[m+1:end])
end
Out[12]:
PowerCurve (generic function with 1 method)
In [13]:
PowerCurve(1000,0.5,600)
Out[13]:
1.3642320780330092e-10
In [14]:
PowerCurve(10,0.5,6)
Out[14]:
0.376953125
In [ ]: