Example 1: Predictivist Bayes Take on a Beta Bernoulli Model.


In [1]:
using Distributions
using Gadfly;

In [2]:
X = [0, 0, 0, 1, 1, 1, 1, 1, 1, 1];

In [3]:
# define the posterior predictive rule for the beta bernoulli in the standard way
PosteriorPredictiveRule1 = function(pastX::Array{Int64})
    α_0 = 1
    β_0 = 1 # define "pseudodata" that will be augmented to any pastX

    if size(pastX,1) > 0
        α_t = α_0 + sum(pastX)
        β_t = β_0 + sum(1 .- pastX)
    else
        α_t = α_0
        β_t = β_0
    end
    
    x_new = rand(Bernoulli(α_t / (α_t + β_t)))
    
    return x_new
end;

In [5]:
# define the terms of the forward simulations
horizon = 1000
nits = 5000;
sims = zeros(Int64,horizon, nits);

In [6]:
# simulate out to the forecast horizon for a number of separate Monte Carlo iterations
for jj = 1:nits
    for ii = 1:horizon
        if ii > 1
            sims[ii,jj] = PosteriorPredictiveRule1(vcat(X,sims[1:(ii-1), jj]))
        else
            sims[ii,jj] = PosteriorPredictiveRule1(X)
        end
    end
end;

In [7]:
# Convert the Bernoulli draws into running means
cumulativesims = cumsum(sims,1)
runningprobs = zeros(Float64, horizon, nits)
sumX = sum(X)
lengthX = size(X,1)
for jj = 1:nits
    for ii = 1:horizon
        runningprobs[ii,jj] = (cumulativesims[ii,jj] + sumX) / (ii + lengthX)
    end
end;

In [8]:
# plot to demonstrate the equivalence of these posterior draws and the analytical Beta(8,4) distribution (Fig 2)
plot(layer(x=runningprobs[horizon,:], Geom.histogram(density=true, bincount=50),order=1),
layer(x=0:.01:1,y=pdf(Beta(8,4),0:0.01:1), Geom.line, Theme(default_color=colorant"grey"),order=2),
Guide.xlabel("θ"), Guide.ylabel("density"))


Out[8]:
θ -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 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 -5 0 5 10 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 density

In [9]:
# plot to demonstrate the simulated sample paths for simulations starting with the observed X (black lines of fig 3)
samplepaths = Layer[]
for ii = 1:20
    push!(samplepaths, layer(x=1:horizon, y=runningprobs[:,ii], Geom.line)[1])
end
plot(samplepaths, Guide.xlabel("iteration"), Guide.ylabel("sample mean"), Coord.Cartesian(ymin=0, ymax=1))


Out[9]:
iteration -1500 -1000 -500 0 500 1000 1500 2000 2500 -1000 -950 -900 -850 -800 -750 -700 -650 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500 1550 1600 1650 1700 1750 1800 1850 1900 1950 2000 -1000 0 1000 2000 -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 -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 sample mean

End of example the first.