"The probability of an event, given a sample space of equiprobable outcomes."
P(event, space)=length(intersect(event,space))//length(space);

D = [1, 2, 3, 4, 5, 6];

even = [2, 4, 6]

P(even, D)

"""The probability of an event, given a sample space of equiprobable outcomes.
    event can be either a set of outcomes, or a predicate (true for outcomes in the event)."""
function P(event, space)
    event_ = such_that(event, space)
    length(intersect(event_,space))//length(space)
end

 #Making use of Julia's multiple dispatch
"The subset of elements in the collection for which the predicate is true."
function such_that(predicate::Function, collection)
    filter(predicate,collection)
end;
"Default return for a collection"
function such_that(event, collection)
    event
end;

even_p(n) = (n % 2 == 0)

such_that(even_p, D)

P(even_p, D)

D12 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]

such_that(even_p, D12)

P(even_p, D12)

S = ["BG", "BB", "GB", "GG"];

two_boys(outcome)=  length(matchall(r"B", outcome)) == 2

older_is_a_boy(outcome)= startswith(outcome,"B");

P(two_boys, such_that(older_is_a_boy, S))

at_least_one_boy(outcome)= 'B' in outcome;

P(two_boys, such_that(at_least_one_boy, S))

such_that(at_least_one_boy, S)

S2b = ["BB/b?", "BB/?b", 
       "BG/b?", "BG/?g", 
       "GB/g?", "GB/?b", 
       "GG/g?", "GG/?g"];

observed_boy(outcome)= 'b' in outcome

such_that(observed_boy, S2b)

P(two_boys, such_that(observed_boy, S2b))

sexesdays = ["$sex$day" 
             for sex in "GB", 
             day in "1234567"]

S3        = ["$older$younger" 
             for older in sexesdays, 
             younger in sexesdays];

@assert length(S3) == (2*7)^2 == 196

print(sort(vec(S3)))

P(at_least_one_boy, S3)

P(at_least_one_boy, S)

P(two_boys, S3)

P(two_boys, S)

P(two_boys, such_that(at_least_one_boy, S3))

P(two_boys, such_that(at_least_one_boy, S))

at_least_one_boy_tues(outcome)= contains(outcome,"B3");

P(two_boys, such_that(at_least_one_boy_tues, S3))

observed_boy_tues(outcome)= contains(outcome,"b3")  

S3b=[]
for children in S3
    for i=1:2
        first_observed=(i==1) ? lowercase(children)[1:2]*"??":"??"*lowercase(children)[3:4]
        push!(S3b,"$children/$first_observed")
    end
end

P(two_boys, such_that(observed_boy_tues, S3b))

"""Display sample space in a table, color-coded: green if event and condition is true; 
   yellow if only condition is true; white otherwise."""
function table(space, n=1, event=two_boys, condition=older_is_a_boy)
    # n is the number of characters that make up the older child.
    olders = sort(unique([outcome[1:n] for outcome in space]))
    html= string("<table>",
                 join([row(older, space, event, condition) for older in olders]),
                "</table>", 
                P(event, such_that(condition, space)))
    display("text/html",html)
end

"Display a row where an older child is paired with each of the possible younger children."
function row(older, space, event, condition)
    thisrow = sort(filter((x)->startswith(x,older),space))
    string("<tr>", join([cell(outcome, event, condition) for outcome in thisrow]),"</tr>")
end

"Display outcome in appropriate color."
function cell(outcome, event, condition)
    color = (event(outcome) && condition(outcome))? "lightgreen": condition(outcome)? "yellow":"ghostwhite"
    return "<td style=\"background-color: $color\">$outcome</td>"      
end;

# Problem 1
table(S, 1, two_boys, older_is_a_boy)

# Problem 2
table(S, 1, two_boys, at_least_one_boy)

# Problem 2
table(S3, 2, two_boys, at_least_one_boy)

# Problem 3
table(S3, 2, two_boys, at_least_one_boy_tues)

B = ["heads/Monday/interviewed", "heads/Tuesday/sleep",
     "tails/Monday/interviewed", "tails/Tuesday/interviewed"];

"Return a predicate that is true of all outcomes that have 'property' as a substring."
T(property)=(outcome)-> contains(outcome,property);

heads = T("heads")
interviewed = T("interviewed")

P(heads, such_that(interviewed, B))

P(heads, B) 

M = ["Car1<Pick1/Open2", "Car1>Pick1/Open3",
     "Car2<Pick1/Open3", "Car2>Pick1/Open3",
     "Car3<Pick1/Open2", "Car3>Pick1/Open2"];

P(T("Car1"), such_that(T("Open3"), M))

P(T("Car2"), such_that(T("Open3"), M))

M2 = ["Car1/Pick1/Open2/Goat", "Car1/Pick1/Open3/Goat",
      "Car2/Pick1/Open2/Car",  "Car2/Pick1/Open3/Goat",
      "Car3/Pick1/Open2/Goat", "Car3/Pick1/Open3/Car"];

P(T("Car1"), such_that(T("Open3/Goat"), M2))

P(T("Car2"), such_that(T("Open3/Goat"), M2))

ProbDist=Dict

"Probability Distribution"
probdist(;entries...)= return normalize(Dict(entries))

"Given a distribution dict, return a version where the values are normalized to sum to 1."
function normalize(dist)
    total = sum(values(dist))
    return Dict(string(e) => dist[e] / total 
        for e in keys(dist))
end;

DK = probdist(GG=121801, GB=126840,
              BG=127123, BB=135138)

"""The probability of an event, given a sample space of equiprobable outcomes. 
    event: a collection of outcomes, or a predicate that is true of outcomes in the event. 
    space: a set of outcomes or a probability distribution of {outcome: frequency} pairs."""
function P(event::Function, space::ProbDist)
    event_ = such_that(event, space)
    return sum([space[v] for v in collect(filter(e-> event(e), keys(space)))])
end
function P(event, space::ProbDist)
    return sum([space[v] for v in collect(filter(e-> e in event, keys(space)))])
end

"""The elements in the space for which the predicate is true.
   If space is a set, return a subset {element,...};
   if space is a dict, return a sub-dict of {element: frequency,...} pairs;
   in both cases only with elements where predicate(element) is true."""
function such_that(predicate::Function, space::ProbDist)
    return normalize(Dict(e=>space[e] for e in filter(predicate,keys(space))))
end;

# Problem 1 in S
P(two_boys, such_that(older_is_a_boy, S))

# Problem 2 in S
P(two_boys, such_that(at_least_one_boy, S))

# Problem 1 in DK
P(two_boys, such_that(older_is_a_boy, DK))

# Problem 2 in DK
P(two_boys, such_that(at_least_one_boy, DK))

"""The joint distribution of two independent probability distributions. 
    Result is all entries of the form {a+b: P(a)*P(b)}"""
joint(A, B)=Dict(a * b => A[a] * B[b] for a in keys(A),b in keys(B));

sexes = probdist(B=51.5, G=48.5)   # Probability distribution over sexes
days  = probdist(L=1, N=4*365)     # Probability distribution over Leap days and Non-leap days
child = joint(sexes, days)        # Probability distribution for one child family
S4    = joint(child, child);     # Probability distribution for two-child family

child

S4

# Problem 4

boy_on_leap_day = T("BL")

P(two_boys, such_that(boy_on_leap_day, S4))

"""Simulate this sequence of events:
- The host randomly chooses a door for the 'car'
- The contestant randomly makes a 'pick' of one of the doors
- The host randomly selects a valid door to be 'opened.' 
- If 'switch' is True, contestant changes 'pick' to the other door
Return true if the pick is the door with the car."""
function monty(switch=true)
    doors  = [1, 2, 3]
    car    = rand(doors)
    pick   = rand(doors)
    opened = rand(filter(d->d != car && d != pick,doors))
    pick = switch? filter(d->d!=pick && d!=opened,doors)[1]: pick 
    pick==car
end;

count(_->monty(true),1:100000)/100000

count(_->monty(false),1:100000)/100000

# The board: a list of the names of the 40 squares
board = """GO   A1 CC1 A2  T1 R1 B1  CH1 B2 B3
           JAIL C1 U1  C2  C3 R2 D1  CC2 D2 D3 
           FP   E1 CH2 E2  E3 R3 F1  F2  U2 F3 
           G2J  G1 G2  CC3 G3 R4 CH3 H1  T2 H2""" |> split

# Lists of 16 community chest and 16 chance cards. See do_card.
CC = append!(["GO", "JAIL"], repeat(["?"],outer=[14]))

CH = append!("GO JAIL C1 E3 H2 R1 R R U -3" |> split,  repeat(["?"],outer=[6]))

"""Simulate given number of steps of monopoly game, 
yielding the name of the current square after each step."""
function monopoly(steps)
    global here
    here = 1
    CC_deck = shuffle(CC)
    CH_deck = shuffle(CH)
    doubles = 0
    function monopolyTask()
        for _=1:steps
            d1, d2 = rand(1:6), rand(1:6)
            goto(here + d1 + d2)
            doubles = (d1 == d2) ? (doubles + 1): 0
            if doubles == 3 || board[here] == "G2J"  
                goto("JAIL")
            elseif startswith(board[here],"CC")
                do_card(CC_deck)
            elseif startswith(board[here],"CH")
                do_card(CH_deck)
            end
            produce(board[here])
        end
    end
    Task(monopolyTask)
end

"Go to destination square (a square number). Update 'here'."
function goto(square::Int)
    global here
    here = (square-1) % length(board)+1
end

"Go to destination square (a square name). Update 'here'."
function goto(square::AbstractString)
    global here
    here = findfirst(board,square)
end

"Take the top card from deck and do what it says."
function do_card(deck)
    global here
    card = pop!(deck)           # The top card
    unshift!(deck,card)           # Move top card to bottom of deck
    if card == "R"|| card == "U" 
        while !startswith(board[here],card)
            goto(here + 1)   # Advance to next railroad or utility
        end
    elseif card == "-3"
        goto(here - 3)       # Go back 3 spaces
    elseif card != "?"
        goto(card)           # Go to destination named on card
    end
end;     

results = collect(monopoly(10^6));

using PyPlot
figure("hist",figsize=(5,3))
ax=axes()
axis([1 ,40, 0 ,70000])
ax[:hist]([findfirst(board,name) for name in results], bins=40);

# We have to implement most_common here at it is not found Julia libraries 
using StatsBase
function most_common(r)
    commonList=collect(zip(board,counts([findfirst(board,name)::Int for name in r])))
    sort!(commonList,by=x->x[2],rev=true)
end
most_common(results)

type Monopoly
    board 
    CC 
    CH 
    here
    function Monopoly()
        this=new()
        this.board= """GO   A1 CC1 A2  T1 R1 B1  CH1 B2 B3
           JAIL C1 U1  C2  C3 R2 D1  CC2 D2 D3 
           FP   E1 CH2 E2  E3 R3 F1  F2  U2 F3 
           G2J  G1 G2  CC3 G3 R4 CH3 H1  T2 H2""" |> split
        this.CC = append!(["GO", "JAIL"], repeat(["?"],outer=[14]))
        this.CH = append!("GO JAIL C1 E3 H2 R1 R R U -3" |> split,  repeat(["?"],outer=[6]))
        shuffle!(this.CC)
        shuffle!(this.CH)
        this.here=1
        return this
    end
end

"""Simulate given number of steps of monopoly game, incrementing counter 
for current square after each step. Return a list of (square, count) pairs in order."""
function simulate(m::Monopoly,steps)
    counter=zeros(Int,40)
     doubles = 0
     for _=1:steps
        d1, d2 = rand(1:6), rand(1:6)
        goto(m, m.here + d1 + d2)
        doubles = (d1 == d2) ? (doubles + 1): 0
        if doubles == 3 || m.board[m.here] == "G2J"  
            goto(m, "JAIL")
        elseif startswith(m.board[m.here],"CC") || startswith(m.board[m.here],"CH")
            do_card(m)
        end
        counter[m.here]+=1
    end
    commonList=collect(zip(m.board,counter))
    sort!(commonList,by=x->x[2],rev=true)
end

"Go to destination square (a square number). Update 'here'."
function goto(m::Monopoly, square::Int)
    m.here = (square-1) % length(m.board)+1
end

"Go to destination square (a square name). Update 'here'."
function goto(m::Monopoly,square::AbstractString)
    m.here = findfirst(m.board,square)
end

"Take the top card from deck and do what it says."
function do_card(m::Monopoly)
    deck= startswith(m.board[m.here],"CC") ? m.CC: m.CH #Which deck based on location
    card = pop!(deck)           # The top card
    unshift!(deck,card)           # Move top card to bottom of deck
    if card == "R"|| card == "U" 
        while !startswith(m.board[m.here],card)
            goto(m,m.here + 1)   # Advance to next railroad or utility
        end
    elseif card == "-3"
        before=m.here
        goto(m,m.here - 3)       # Go back 3 spaces
        if m.here==0
            println(before, " ",card)
        end
    elseif card != "?"
        goto(m,card)           # Go to destination named on card
    end
 
end
                    
simulate(Monopoly(),10^6)

"Return the probability distribution for the St. Petersburg Paradox with a limited bank."
function st_pete(limit)
    P = Dict()     # The probability distribution
    pot = 2    # Amount of money in the pot
    pr  = 1/2 # Probability that you end up with the amount in pot
    while pot < limit 
        P[pot] = pr
        pot, pr = pot * 2, pr / 2
    end
    P[limit] = pr * 2               # pr * 2 because you get limit for heads or tails
    assert(sum(values(P)) == 1.0)
    sort(collect(zip(keys(P),values(P))))
end

StP = st_pete(10^9)

"The expected value of a probability distribution."
EV(P) =sum([v[1] * v[2] for v in P]);

EV(StP)

"The value of money: only half as valuable after you already have enough."
function util(dollars, enough=1000)
    if dollars < enough
        return dollars
    else
        return enough + util((dollars-enough)/2., enough*2)
    end
end;

for d=2:10
    println(lpad(10^d,15)," \$ = ",lpad(round(Int,util(10^d)),10), " util")
end

figure("Value of Money",figsize=(5,3))
plot([util(x) for x=1000:1000:10000000])
println("Y axis is util(x); x axis is in thousands of dollars.")

"The expected utility of a probability distribution, given a utility function."
EU(P, U)= sum([e[2] * U(e[1]) for e in P]);

EU(StP, util)

flip()=rand(["head", "tail"])

"Simulate one round of the St. Petersburg game, and return the payoff."
function simulate_st_pete(limit=10^9)
    pot = 2
    while flip() == "head"
        pot = pot * 2
        if pot > limit
            return limit
        end
    end
    return pot
end;

srand(1234)
# With Julia to can take a few million samples at get the results instantly 
results = sort([(z[1],z[2]) for z in proportionmap([simulate_st_pete() for _=1:100_000])], by=x->x[1])

EU(results, util) 

EV(results)

"For each element in the iterable, yield the mean of all elements seen so far." 
function running_averages(iterable)
    total, n = 0, 0
    function run_avg_task()
        for x in iterable
            total, n = total + x, n + 1
            produce(total / n)
            end
    end
    Task(run_avg_task)
end

"Plot the running average of calling the function n times."
function plot_running_averages(fn, n)
    plot(collect(running_averages([fn() for _=1:n])))
end;

srand(555)
figure("Running Average",figsize=(5,3))
for i=1:10
    plot_running_averages(simulate_st_pete, 100000)
end
axis([1 ,100_000, 2, 140]);