Monte Carlo weather forecasting competition

forecastcompetition.jl

using Plots
using Distributions

"""Continuous ground truth."""
function generate_ground_truth(rng::Tuple{Integer, Integer}, n=14)
    (rand(Int, n) .% (rng[2]-rng[1]+1) .+ rng[1])
end

"""Bernoulli ground truth."""
function generate_ground_truth(prob1::Real, n=14)
    rand(Bernoulli(prob1), n)
end

"""Continuous predictions."""
function generate_forecasts(groundtruth::Vector, σ=3., m=5)
    hcat([randn(m).*σ .+ gt for gt in groundtruth]...)
end

"""Bernoulli predictions."""
function generate_forecasts_bn(groundtruth::Vector{T}, p, m=5) where {T <: Integer}
    hcat([rand(Bernoulli(p), m) for gt in groundtruth]...)
end

linerr(x) = abs(x)
sqerr(x) = x^2
normpdf(x, σ) = pdf(Normal(0, σ), x)
bernoulliscore(x, factor=13) = factor*abs(x)

function judge(groundtruth::Vector, forecasts::Matrix, score=normpdf)
    if length(groundtruth) != size(forecasts, 2)
        error("judge(): need as many forecasts as groundtruths")
    end
    sum(score.(forecasts .- groundtruth'); dims=2)
end

function simulate(days=14, participants=5)
    n = days
    m = participants

    scores = zeros(m)

    # Temp low
    gtl = generate_ground_truth((40, 60), n)
    fclow = generate_forecasts(gtl, 3., m)
    scores .+= judge(gtl, fclow, sqerr)#x -> 1/normpdf(x, 3.))

    p1 = plot(legend=false, xlabel="day", ylabel="temp")
    bar!(p1, gtl)
    scatter!(fclow')

    # Temp high
    gth = generate_ground_truth((70, 100), n)
    fch = generate_forecasts(gth, 3., m)
    scores .+= judge(gth, fch, sqerr) #x -> 1/normpdf(x, 3.))

    # Gusts
    gtg = generate_ground_truth(.2, n)
    fcgust = generate_forecasts_bn(gtg, .3, m)
    scores .+= judge(gtg, fcgust, bernoulliscore)

    scores, [p1]
end

function plot_many(days, participants, experiments)
    runs = [(simulate(days, participants)) for i in 1:experiments]
    plots = [r[2] for r in runs]
    runs = sort.([r[1] for r in runs])

    p = plot(legend=false, yscale=:log10,
        xlabel="participant (sorted by rank, lower is better)",
        ylabel="log10(score) (lower is better)")
    plot!(p, runs)

    vcat(p, plots...)
end