Sampling#

What does in mean to sample from a distribution?

It is useful to consider two separate cases: a discrete distribution and a continuous one. Lets start with the discrete distribution.

mysample()=if rand()<0.5 ; -1; else ;1; end

mysample (generic function with 1 method)

using Plots
mysamples=[mysample() for j in 1:100]
histogram(mysamples)

using StatsBase

sample([1,-1],Weights([0.5,0.5]))

-1

Weights([0.5,0.5])

2-element Weights{Float64, Float64, Vector{Float64}}:
 0.5
 0.5

rand()

0.04766349983427909

rand(5)

5-element Vector{Float64}:
 0.22263402030579593
 0.830845210182813
 0.061544518835685236
 0.6115764528499936
 0.5453439465782544

rand(2,3)

2×3 Matrix{Float64}:
 0.107404  0.429609  0.996006
 0.25372   0.181249  0.207257

2*rand()-1 # gives random number from -1 to 1

-0.6613577730458822

ylist=rand(1000)
xlist=sqrt.(ylist)
histogram(xlist)

rads=rand(1000)
angles=2*pi*rand(1000)
pnts=[(r*cos(ϕ),r*sin(ϕ)) for (r,ϕ) in zip(rads,angles)]
scatter(pnts,aspect_ratio=:equal)

sqrads=rand(1000)
angles=2*pi*rand(1000)
rads=sqrt.(sqrads)
pnts=[(r*cos(ϕ),r*sin(ϕ)) for (r,ϕ) in zip(rads,angles)]
scatter(pnts,aspect_ratio=:equal)

function boxmuller(n)
    halfn=Int64(ceil(n/2))
    svals=rand(halfn)
    phivals=2*pi*rand(halfn)
    scaledsvals= @. sqrt(-log(svals))
    vcat(scaledsvals.*cos.(phivals),scaledsvals.*cos.(phivals))[1:n]
end

boxmuller (generic function with 1 method)

histogram(boxmuller(50000),normed=true,label="")
plot!(x->exp(-x^2)/sqrt(pi),-2,2,label="",xlabel="x",ylabel="p")

using Distributions

d=Normal()

Normal{Float64}(μ=0.0, σ=1.0)

histogram(rand(d,50000),label="",normed=true)
plot!(x->exp(-x^2/2)/sqrt(2*pi),-2,2,label="",xlabel="x",ylabel="p")

Rejection sampling#

Another way to sample from a non-uniform distribution is to use rejection sampling. Lets return to our example of

What we will do is uniformly sample $x$ from $[0,1]$. Then we will roll a random number. If that random number is less than $x$, we keep it, otherwise we start over and draw a new $x$.

function rsample()
    x=0.
    while true
        x=rand()
        y=rand()
        if y<x
            break
        end
    end
    return x
end

rsample(n)=[rsample() for j in 1:n]

rsample (generic function with 2 methods)

histogram(rsample(10000),normed=true)

The reason this works is we can write

The probability of selecting $x$ is uniform, while the acceptance probability is proportional to $x$ – which gives the desired result.

Markov Chains#

So far all of our sampling methods have given us independent samples. There are also approaches to finding series of samples, where the individual realizations are correlated, but after a large number of steps you have an unbiased sample.

A great example is shuffling cards. Every time you do a riffle shuffle, you somewhat mix up the cards, but they are not fully randomized. However, if you do many riffles, you end up with a truly random set of cards.

A classic example for thinking about this is the Ehrenfest urn model, also known as the dog and flea model. Suppose you have two dogs, which have in total $N$ fleas on them. At each time step one of the fleas hops from one dog to the other.

Suppose there are $j$ fleas on the left dog, and $N-j$ on the right dog. You would expect that the probability of a hop from left to right would be proportional to the number of fleas on the left dog. Similarly the probability of a hop from right to left is proportional to the number on the right one. We will refer to the configuration of fleas as a state and label them by the number of fleas on the left dog.

We say that in each step the probability of going from state $j$ to $j+1$ is ${\mathrm{\Gamma }}_{j\to j+1}=j/N$ and the probability of going from $j$ to $j-1$ is ${\mathrm{\Gamma }}_{j\to j-1}=(N-j)/N$. You can verify that these probabilities add to 1. Lets program this model.

function simulate_urn(initialj,N,numsteps)
    j=initialj
    path=[j]
    for i in 1:numsteps
        gammaplus=(N-j)/N
        if rand()<gammaplus
            j+=1
        else
            j-=1
        end
        push!(path,j)
    end
    path
end   

simulate_urn (generic function with 1 method)

Lets now simulate this for 100 steps, repeating the process 1000 times.

plt=plot(xlabel="t",ylabel="j")
for j in 1:50
    plot!(simulate_urn(50,100,1000),label="")
end
display(plt)

As you can see, things settle into some sort of steady state. Here is the distribution at the final time

data=hcat([simulate_urn(50,100,1000) for j in 1:2000]...)
histogram(data[end,:],bins=-0.5:1:101.5,normed=true)

We can get better statistics if we bin all of the times:

histogram(vcat(data[100:end,:]...),bins=-0.5:1:100.5,normed=true)

How can we understand this distribution?

Let $p_j(t)$ be the probability of site $j$ at time $t$. By our rules,

It is useful to write this in terms of the change in probability:

The terms in square brackets have a physical meaning: They are the probability fluxes. They are the net flow of probability between two sites.

After long time we reach a steady state, and these fluxes should vanish:

In our particular case ${\mathrm{\Gamma }}_{j\to j+1}=j/N$ and ${\mathrm{\Gamma }}_{j+1\to j}=(N-j-1)/N$ and hence

One can show that the solution to this recurrance relationship is the binomial distribution

which we encountered earlier with our coin flip example. Thus we can use this stochastic sequence to sample from the binomial distribution.

histogram(vcat(data[100:end,:]...),bins=-0.5:1:100.5,normed=true)
jlist=0:100
plist=[binomial(big(100),big(j))/(2.)^100 for j in jlist]
plot!(jlist,plist,label="")

Here I used big to convert integers into BigInt objects – which are a form of integer which can be arbitrarily large

factorial(big(100))

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

factorial(100)

OverflowError: 100 is too large to look up in the table; consider using `factorial(big(100))` instead

Stacktrace:
 [1] factorial_lookup
   @ ./combinatorics.jl:19 [inlined]
 [2] factorial(n::Int64)
   @ Base ./combinatorics.jl:27
 [3] top-level scope
   @ In[26]:1

The Ehrenfest model is an example of a Markov chain. A Markov chain is the simplest form of stochastic process. It consists of a set of states $i$ (discrete in this case, but they could be continuous), and a set of transition rates ${\mathrm{\Gamma }}_{i\to j}$. If at time $t$ you are in state $i$, then you transition to state $j$ with probability ${\mathrm{\Gamma }}_{i\to j}$. Since probability must be normalized

though I should note that you are generically allowed to stay put: ${\mathrm{\Gamma }}_{i\to i}$ can be non-zero.

Typically a Markov chain will reach a steady-state distribution. By definition, the steady state distribution obeys a ballance condition

The net flux in to a state is equal to the net flux out. In most cases that steady state distribution satisfies a stronger condition – detailed balance – meaning that the net flux between any two state vanishes:

Using Markov chains to sample arbitrary distributions#

An important application of Markov chains is we can use them to sample arbitrary distributions. For example, suppose we want to sample $p_j\propto j$ for $j=1,2,\cdots 100$. We could set up a Markov process where ${\mathrm{\Gamma }}_{j\to j+1}/{\mathrm{\Gamma }}_{j+1\to j}=(j+1)/j$.

One algorithm for coming up with the transition rates is the Metropolis algorithm. It explicitly constructs ${\mathrm{\Gamma }}_{i\to j}$ for $i\ne j$, then implicitly ${\mathrm{\Gamma }}_{i\to i}=1-\sum _j{\mathrm{\Gamma }}_{i\to j}$.

We break ${\mathrm{\Gamma }}_{i\to j}$ into two pieces: a proposal ${\mathrm{\Pi }}_{i\to j}$ and an acceptance $A_{i\to j}$. Given our current state $i$, we choose a random move to $j$ with probability distribution ${\mathrm{\Pi }}_{i\to j}$, then accept that move with probability $A_{i\to j}$. If we don’t accept, then we stay put. The proposals are a probability distribution, so we need $\sum _j{\mathrm{\Pi }}_{i\to j}=1$.

In this case, lets keep the proposal simple. We will only propose moves which change $j$ by 1. Lets take ${\mathrm{\Pi }}_{j\to j+1}={\mathrm{\Pi }}_{j\to j-1}=1/2$. To maintain normalization, on the edges we will also need to propose “staying put” as a move: ${\mathrm{\Pi }}_{1\to 1}={\mathrm{\Pi }}_{100\to 100}=1/2$. Thus every proposed move will be equally likely.

Next, we choose $A_{i\to j}$ to satisfy detailed balance. One way to do this is take

In our case $A_{j\to j+1}=1$, and $A_{j+1\to j}=j/(j+1)$.

Lets code that up, passing the (non-normalized) distribution to the function

function markov_walk(initialj,distribution,numsteps)
    j=initialj
    path=[j]
    for i in 1:numsteps
        # randomly propose left or right
        if rand()<1/2 # propose left
            if j==1 || rand()>min(1,distribution[j-1]/distribution[j])
                push!(path,j)
                continue # go to next step of i loop
            end
            j-=1
            push!(path,j)
            continue
        end
        # We have chosen right
        if j==length(distribution) || 
                rand()>min(1,distribution[j+1]/distribution[j])
            push!(path,j)
            continue
        end
        j+=1
        push!(path,j)
    end
    path
end   

markov_walk (generic function with 1 method)

w=markov_walk(1,collect(1:100),1000000)
plot(w[1:1000])

plot(w[1:100:end])

histogram(w[1000:100:end],bins=-0.5:1:100.5)