diff --git a/Ch6.ipynb b/Ch6.ipynb deleted file mode 100644 index 55f06ec..0000000 --- a/Ch6.ipynb +++ /dev/null @@ -1,1177 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Setup" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# parallel computing\n", - "using Distributed\n", - "addprocs(12);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# needed packages\n", - "@everywhere using Distributions, Bootstrap, Statistics, LinearAlgebra, SharedArrays\n", - "using DataFrames, HypothesisTests, DelimitedFiles, StatsBase, Colors, Gadfly" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# default graphics\n", - "Gadfly.push_theme(:default)\n", - "set_default_plot_size(9inch, 9inch/MathConstants.golden)\n", - "\n", - "function gen_brew_colors(n) # to create your own colors, here based on one of the brewer series\n", - " cs = distinguishable_colors(n, \n", - " [colorant\"#66c2a5\", colorant\"#fc8d62\", colorant\"#8da0cb\", colorant\"#e78ac3\",\n", - " colorant\"#a6d854\", colorant\"#ffd92f\", colorant\"#e5c494\", colorant\"#b3b3b3\"],\n", - " lchoices=Float64[58, 45, 72.5, 90],\n", - " transform=c->deuteranopic(c, 0.1),\n", - " cchoices=Float64[20,40],\n", - " hchoices=[75,51,35,120,180,210,270,310]\n", - " )\n", - " convert(Vector{Color}, cs)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# set parameters, define priors, etc.\n", - "@everywhere begin \n", - " const numb_hyp = 11\n", - " const numb_toss = 500\n", - " const numb_sim = 1000\n", - " const prior = fill(Float32(1/numb_hyp), numb_hyp)\n", - " const likelihood_heads = range(0f0, stop=1, length=numb_hyp)\n", - " const likelihood_tails = range(1f0, stop=0, length=numb_hyp)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere datFunc(bias) = rand(Bernoulli(bias), numb_toss)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Update rules" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Bayes' rule\n", - "@everywhere function b_upd(probs::Array{Float32,1}, dat::Array{Bool,1}, toss_num::Int64)\n", - " if dat[toss_num] == true\n", - " @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)\n", - " else\n", - " @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# EXPL\n", - "@everywhere function expl_upd(probs::Array{Float32,1}, dat::Array{Bool, 1}, toss_num::Int64, bonus::Float32=0.1)\n", - " val::Float32 = mean(dat[1:toss_num]) * 10 + 1\n", - " vec::Array{Float32,1} = if dat[toss_num] == true\n", - " @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)\n", - " else\n", - " @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)\n", - " end\n", - "\n", - " if val % 1 == .5\n", - " vec[floor(Int64, val)] += .5*bonus\n", - " vec[ceil(Int64, val)] += .5*bonus\n", - " else\n", - " vec[round(Int64, val, RoundNearestTiesAway)] += bonus\n", - " end\n", - "\n", - " return vec / (1 + bonus)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Good's rule\n", - "@everywhere function good_bonus(probs::Array{Float32,1}, res::Bool, λ=2) # with λ=2, we obtain the rule L2 from Douven and Schupbach, Frontiers ...\n", - "\n", - " pE::Float32 = res == true ? dot(probs, likelihood_heads) : dot(probs, likelihood_tails)\n", - " gb::Array{Float32,1} = res == true ? log.(likelihood_heads ./ pE) : log.(likelihood_tails ./ pE)\n", - "\n", - " function rsc(i)\n", - " if i >= 0\n", - " 1 - exp(2λ^2 * -i^2)\n", - " else\n", - " -1 + exp(2λ^2 * -i^2)\n", - " end\n", - " end\n", - "\n", - " return map(rsc, gb)\n", - "\n", - "end\n", - "\n", - "@everywhere function good_upd(probs::Array{Float32,1}, dat::Array{Bool,1}, toss_num::Int64, γ::Float32=0.1)\n", - "\n", - " res::Bool = dat[toss_num]\n", - "\n", - " probvec::Array{Float32,1} = if res == true\n", - " @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)\n", - " else\n", - " @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)\n", - " end\n", - "\n", - " goodvec::Array{Float32,1} = probvec + γ .* (probvec .* good_bonus(probs, res))\n", - "\n", - " return goodvec / sum(goodvec)\n", - "\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Popper's rule\n", - "@everywhere function pop_bonus(probs::Array{Float32,1}, res::Bool)\n", - "\n", - " pE::Float32 = res == true ? dot(probs, likelihood_heads) : dot(probs, likelihood_tails)\n", - " pb::Array{Float32, 1} = res == true ? (likelihood_heads .- pE) ./ (likelihood_heads .+ pE) : (likelihood_tails .- pE) ./ (likelihood_tails .+ pE)\n", - "\n", - " end\n", - "\n", - "@everywhere function pop_upd(probs::Array{Float32,1}, dat::Array{Bool, 1}, toss_num::Int64, γ::Float32=0.1)\n", - "\n", - " res::Bool = dat[toss_num]\n", - "\n", - " probvec::Array{Float32,1} = if res == true\n", - " @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)\n", - " else\n", - " @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)\n", - " end\n", - "\n", - " popvec::Array{Float32,1} = probvec + γ .* (probvec .* pop_bonus(probs, res))\n", - "\n", - " return popvec / sum(popvec)\n", - "\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Survival distributions" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Weibull distribution" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "weibDistr = plot(\n", - " [x->cdf(Weibull(1, i), x) for i in 50:50:250][:],\n", - " color=[\"Weibull(1, $i)\" for i in 50:50:250][:],\n", - " 0, 100,\n", - " Guide.colorkey(title=\"Distribution\"),\n", - " Guide.xlabel(\"Time\"),\n", - " Guide.ylabel(\"Probability of death\"),\n", - " Guide.title(\"CDFs of Weibull distributions\"),\n", - " Scale.color_discrete_manual(gen_brew_colors(5)...),\n", - " style(line_width=2pt, minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "weibInter = plot(\n", - " [x->(1 + cdf(Weibull(1, 50), x))/2,\n", - " x->cdf(Weibull(1, 50), x),\n", - " x->cdf(Weibull(1, 50), x)/2,\n", - " ],\n", - " color=[\"wrong\", \"none\", \"right\"],\n", - " 0, 100,\n", - " Guide.colorkey(title=\"Intervention\"),\n", - " Guide.xlabel(\"Time\"),\n", - " Guide.ylabel(\"Probability of death\"),\n", - " Guide.title(\"Effect of intervention: Weibull(1, 50)\"),\n", - " Scale.color_discrete_manual(gen_brew_colors(3)...),\n", - " style(line_width=2pt, minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Gamma distribution" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "gamDistr = plot(\n", - " [x->cdf(Gamma(10, i), x) for i in [10, 12, 14, 16]][:],\n", - " color=[\"Γ(10, $i)\" for i in [10, 12, 14, 16]][:],\n", - " 0, 500,\n", - " Guide.colorkey(title=\"Distribution\"),\n", - " Guide.xlabel(\"Time\"),\n", - " Guide.ylabel(\"Probability of death\"),\n", - " Guide.title(\"CDFs of Gamma distributions\"),\n", - " Scale.color_discrete_manual(gen_brew_colors(4)...),\n", - " style(line_width=2pt, minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "gammaInter = plot(\n", - " [x->(1 + cdf(Gamma(10, 16), x))/2,\n", - " x->cdf(Gamma(10, 16), x),\n", - " x->cdf(Gamma(10, 16), x)/2,\n", - " ],\n", - " color=[\"wrong\", \"none\", \"right\"],\n", - " 0, 500,\n", - " Guide.colorkey(title=\"Intervention\"),\n", - " Guide.xlabel(\"Time\"),\n", - " Guide.ylabel(\"Probability of death\"),\n", - " Guide.title(\"Effect of intervention for Γ(10, 16)\"),\n", - " Scale.color_discrete_manual(gen_brew_colors(3)...),\n", - " style(line_width=2pt, minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Running the simulations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere const numb_agents = 200\n", - "@everywhere const numb_generations = 100" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# starting position: 50 Bayesians, and 50 agents per other group (EXPL, Good's rule, Popper's rule), with varying values for c (varying between 0 and 0.25)\n", - "@everywhere const numb_agents = 200\n", - "groupID = repeat(1.0:4.0, inner=div(numb_agents, 4))\n", - "population_start = vcat(fill(0, div(numb_agents, 4)), rand(Uniform(0, .25), 3*div(numb_agents, 4)))\n", - "pop_start = Array{Float32,2}(hcat(groupID, population_start));" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere function survWei(upds::Array{Float32,2}, # modeling probability of death, based on Weibull distribution\n", - " hyp::Int64,\n", - " a::Float64, \n", - " b::Float64,\n", - " shape::Float64=rand(Uniform(.5, 5)), \n", - " scale::Float64=rand(Uniform(50, 250)), \n", - " thresh::Float64=.9)\n", - " \n", - " t = something(findfirst(upds .> thresh), (numb_toss, 0)) # where in the matrix with probability updates do we find the first value above thresh?\n", - " c = t[2]\n", - " p = t[1]\n", - " \n", - " # cdf(Weibull(shape, scale), p) below gives the probability of death at the relevant time\n", - "\n", - " if c == hyp\n", - " 1 - (cdf(Weibull(shape, scale), p) / a) # probability goes down if right intervention is made (which is made when the truth is assigned a probability above thresh)\n", - " elseif c == 0\n", - " 1 - cdf(Weibull(shape, scale), numb_toss + 1) # if no intervention is made, output survival probability at last time step\n", - " else\n", - " (1 + (b - 1) * cdf(Weibull(shape, scale), p)) / b # probability goes down if wrong intervention is made (which happens if a false hypothesis is assigned a probabilty above thresh)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere function survGam(upds::Array{Float32,2}, # modeling probability of death, based on Gamma distribution\n", - " hyp::Int64,\n", - " a::Float64, \n", - " b::Float64,\n", - " shape::Float64=rand(Uniform(10, 16)), \n", - " scale::Float64=rand(Uniform(10, 16)), \n", - " thresh::Float64=.9)\n", - " \n", - " t = something(findfirst(upds .> thresh), (numb_toss, 0)) # where in the matrix with probability updates do we find the first value above thresh?\n", - " c = t[2]\n", - " p = t[1]\n", - " \n", - " if c == hyp\n", - " 1 - (cdf(Gamma(shape, scale), p) / a) # the probability goes down if the right intervention is made (and the right intervention is made if the truth is assigned a probability above thresh)\n", - " elseif c == 0\n", - " 1 - cdf(Gamma(shape, scale), numb_toss + 1) # if no intervention is made, output survival probability at last time step\n", - " else\n", - " (1 + (b - 1) * cdf(Gamma(shape, scale), p)) / b # the probability goes down if the wrong intervention is made (which happens if a false hypothesis is assigned a probabilty above thresh)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere function patient(rule_index::Float32, c_value::Float32, dist::Function)\n", - "\n", - " rand_hyp::Int64 = rand(1:11) # pick α hypothesis (\"what's wrong with the patient\")\n", - " right = rand(Uniform(1, 10)) # effect of right intervention\n", - " wrong = rand(Uniform(1, 10)) # effect of wrong intervention\n", - "\n", - " data::Array{Bool, 1} = datFunc((rand_hyp - 1) / (numb_hyp - 1)) # generate synthetic data for this pick (the test results for the patient)\n", - " \n", - " updates = Array{Float32,2}(undef, numb_toss + 1, numb_hyp) # initialize array for probabilities\n", - "\n", - " updates[1, :] = prior # set prior\n", - "\n", - " if rule_index == 1.0f0\n", - " @fastmath @inbounds for t in 1:numb_toss # generate updates\n", - " updates[t + 1, :] = b_upd(updates[t, :], data, t)\n", - " end\n", - " elseif rule_index == 2.0f0\n", - " @fastmath @inbounds for t in 1:numb_toss # generate updates\n", - " updates[t + 1, :] = expl_upd(updates[t, :], data, t, c_value)\n", - " end\n", - " elseif rule_index == 3.0f0\n", - " @fastmath @inbounds for t in 1:numb_toss # generate updates\n", - " updates[t + 1, :] = good_upd(updates[t, :], data, t, c_value)\n", - " end\n", - " else\n", - " @fastmath @inbounds for t in 1:numb_toss # generate updates\n", - " updates[t + 1, :] = pop_upd(updates[t, :], data, t, c_value)\n", - " end\n", - " end\n", - "\n", - " return dist(updates, rand_hyp, right, wrong)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#= tests doctor on 100 patients and calculates average survival score obtained by doctor, so\n", - "average probability that patients would survive =#\n", - "@everywhere @inbounds function avScore(rule_index::Float32, c_value::Float32, dist::Function)\n", - " tot = @distributed (+) for i in 1:100\n", - " patient(rule_index, c_value, dist)\n", - " end\n", - " return tot / 100\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#= tests all doctors in a population and selects the 50 percent fittest; outputs those\n", - "as well as a copy of each of them =#\n", - "function population_upd_rep(pop::Array{Float32,2}, dist::Function)\n", - " agent_scores = SharedArray{Float32,1}(numb_agents)\n", - " @sync @distributed for i in 1:numb_agents\n", - " @inbounds agent_scores[i] = avScore(pop[i, :]..., dist)\n", - " end\n", - " best_index = findall(agent_scores .>= Statistics.median(agent_scores))\n", - " best = pop[best_index[1:Int(numb_agents/2)], :]\n", - " return vcat(best, best), agent_scores\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#= runs the foregoing function for 100 generations (together with the starting population,\n", - "we have 101 generations in total), registering for each generation all relevant agent properties,\n", - "so the type the doctor belongs to, the bonus value attributed by the doctor, as well as the doctor's\n", - "fitness score =#\n", - "pop_upd_c_a = Array{Float32,3}(undef, numb_agents, 2, numb_generations + 1)\n", - "pop_upd_f = Array{Float32,2}(undef, numb_agents, numb_generations + 1)\n", - "\n", - "pop_upd_c_a[:, :, 1] = pop_start\n", - "\n", - "@inbounds for i in 1:numb_generations\n", - " pop_upd_c_a[:, :, i + 1], pop_upd_f[:, i] = population_upd_rep(pop_upd_c_a[:, :, i], survWei)\n", - "end\n", - "\n", - "@inbounds for i in 1:numb_agents\n", - " pop_upd_f[i, numb_generations + 1] = avScore((pop_upd_c_a[i, 1:2, numb_generations + 1])..., survWei)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# extracts agent type and bonus values, for easy separate storage\n", - "res_a = Array{Int32,2}(undef, numb_agents, numb_generations + 1)\n", - "res_c = Array{Float32,2}(undef, numb_agents, numb_generations + 1)\n", - "\n", - "for i in 1:numb_generations + 1\n", - " res_a[:, i], res_c[:, i] = pop_upd_c_a[:, 1, i], pop_upd_c_a[:, 2, i]\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "run(`mkdir data`);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# stores the relevant data for all generations\n", - "writedlm(\"data/weib_agent_type1.txt\", res_a)\n", - "writedlm(\"data/weib_c_value1.txt\", res_c)\n", - "writedlm(\"data/weib_fit1.txt\", pop_upd_f)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The above allows running one simulation at a time, the relevant output of which can then be stored. The function below runs 100 simulations and stores the relevant output of each. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function sim_run(dist::Function)\n", - " k = 1\n", - " while k < 101\n", - " groupID = repeat(1.0:4.0, inner=div(numb_agents, 4))\n", - " population_start = vcat(fill(0, div(numb_agents, 4)), rand(Uniform(0, .25), 3*div(numb_agents, 4)))\n", - " pop_start = Array{Float32,2}(hcat(groupID, population_start))\n", - " pop_upd_c_a = Array{Float32,3}(undef, numb_agents, 2, numb_generations + 1)\n", - " pop_upd_f = Array{Float32,2}(undef, numb_agents, numb_generations + 1)\n", - " pop_upd_c_a[:, :, 1] = pop_start\n", - "\n", - " @inbounds for i in 1:numb_generations\n", - " pop_upd_c_a[:, :, i + 1], pop_upd_f[:, i] = population_upd_rep(pop_upd_c_a[:, :, i], dist)\n", - " end\n", - "\n", - " @inbounds for i in 1:numb_agents\n", - " pop_upd_f[i, numb_generations + 1] = avScore(pop_upd_c_a[i, 1:2, numb_generations + 1]..., dist)\n", - " end\n", - "\n", - " res_a = Array{Int32,2}(undef, numb_agents, numb_generations + 1)\n", - " res_c = Array{Float32,2}(undef, numb_agents, numb_generations + 1)\n", - "\n", - " @inbounds for i in 1:(numb_generations + 1)\n", - " res_a[:, i], res_c[:, i] = pop_upd_c_a[:, 1, i], pop_upd_c_a[:, 2, i]\n", - " end\n", - "\n", - " writedlm(\"data/weib_agent_type$k.txt\", res_a)\n", - " writedlm(\"data/weib_c_value$k.txt\", res_c)\n", - " writedlm(\"data/weib_fit$k.txt\", pop_upd_f)\n", - "\n", - " population_start = nothing\n", - " pop_start = nothing\n", - " pop_upd_c_a = nothing\n", - " pop_upd_f = nothing\n", - " res_a = nothing\n", - " res_c = nothing\n", - " GC.gc()\n", - "\n", - " k += 1\n", - " end\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sim_run(survWei)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plot counts of types per generation for one simulation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fullType = readdlm(\"data/weib_agent_type1.txt\");" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ks = [ keys(sort(countmap(fullType[:,i]))) for i in 1:numb_generations + 1 ]\n", - "vls = [ values(sort(countmap(fullType[:,i]))) for i in 1:numb_generations + 1 ]\n", - "\n", - "group = []\n", - "freq = []\n", - "gen = []\n", - "\n", - "for i in 1:101\n", - " append!(group, collect(ks[i]))\n", - " append!(freq, collect(vls[i]))\n", - " append!(gen, fill(i, length(collect(ks[i]))))\n", - "end\n", - "\n", - "bar_df = convert(DataFrame, hcat(group, freq, gen))\n", - "\n", - "bar_df[!, :x4] = map(bar_df[!, :x1]) do x\n", - " if x == 1\n", - " return \"Bayes\"\n", - " elseif x == 2\n", - " return \"EXPL\"\n", - " elseif x == 3\n", - " return \"Good\"\n", - " else\n", - " return \"Popper\"\n", - " end\n", - "end\n", - "\n", - "rename!(bar_df, [Symbol(\"$i\") for i in [\"Group\", \"Count\", \"Generation\", \"Rule\"]]);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(bar_df, x=:Generation, y=:Count, color=:Rule, Geom.bar(position=:stack),\n", - " Coord.cartesian(xmin=1, xmax=numb_generations + 1),\n", - " Guide.colorkey(title=\"Rule\"),\n", - " Scale.color_discrete_manual(gen_brew_colors(4)...),\n", - " style(minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " grid_color=colorant\"#222831\",\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Visualize full results" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "full_results_type = Array{Int32, 3}(undef, 200, 101, 100)\n", - "\n", - "for i in 1:100\n", - " full_results_type[:, :, i] = readdlm(\"data/weib_agent_type$i.txt\")\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "percentages = Array{Float64, 3}(undef, 100, 101, 4)\n", - "\n", - "for k in 1:4\n", - " for j in 1:101\n", - " percentages[:, j, k] = [ length(findall(full_results_type[:, j, i] .== k)) / 200 for i in 1:100 ]\n", - " end\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere function bs(x, n=1000)\n", - " bootstrap(mean, x, BasicSampling(n))\n", - "end\n", - "\n", - "bayesBoot = SharedArray{Float64}(101, 3);\n", - " \n", - "@distributed for i in 1:101\n", - " bayesBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 1]), BasicConfInt(.95))[1]...]\n", - "end\n", - "\n", - "explBoot = SharedArray{Float64}(101, 3);\n", - " \n", - "@distributed for i in 1:101\n", - " explBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 2]), BasicConfInt(.95))[1]...]\n", - "end\n", - "\n", - "goodBoot = SharedArray{Float64}(101, 3);\n", - " \n", - "@distributed for i in 1:101\n", - " goodBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 3]), BasicConfInt(.95))[1]...]\n", - "end\n", - "\n", - "popperBoot = SharedArray{Float64}(101, 3);\n", - " \n", - "@distributed for i in 1:101\n", - " popperBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 4]), BasicConfInt(.95))[1]...]\n", - "end;" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Generation = repeat(collect(1:101), outer=4)\n", - "Rule = repeat([\"Bayes\", \"EXPL\", \"Good\", \"Popper\"], inner=101)\n", - "type_df = convert(DataFrame, hcat(vcat(bayesBoot, explBoot, goodBoot, popperBoot), Generation, Rule))\n", - "rename!(type_df, [Symbol(\"$i\") for i in [\"y\", \"ymin\", \"ymax\", \"Generation\", \"Rule\"]]);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(type_df, x=:Generation, y=:y, ymin=:ymin, ymax=:ymax, color=:Rule, Geom.line, Geom.ribbon,\n", - " Guide.ylabel(\"Average percentage\"),\n", - " Coord.cartesian(xmin=1, xmax=101, ymin=-.001),\n", - " Scale.color_discrete_manual(gen_brew_colors(4)...),\n", - " style(line_width=2pt, lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2),\n", - " minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We do the same for the bonus values." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "full_results_c = Array{Float32, 3}(undef, 200, 101, 100)\n", - "\n", - "for i in 1:100\n", - " full_results_c[:, :, i] = readdlm(\"data/weib_c_value$i.txt\")\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "expl_c_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_c[:, k, i][(full_results_type[:, k, i] .== 2)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " expl_c_res[k, :] = vcat(mean(cr), [confint(OneSampleTTest(cr))...])\n", - "end\n", - "\n", - "good_c_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_c[:, k, i][(full_results_type[:, k, i] .== 3)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " good_c_res[k, :] = vcat(mean(cr), [confint(OneSampleTTest(cr))...])\n", - "end\n", - "\n", - "pop_c_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_c[:, k, i][(full_results_type[:, k, i] .== 4)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " pop_c_res[k, :] = vcat(mean(cr), [confint(OneSampleTTest(cr))...])\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Generation = repeat(collect(1:101), outer=3)\n", - "Rule = repeat([\"EXPL\", \"Good\", \"Popper\"], inner=101)\n", - "cval_df = convert(DataFrame, hcat(vcat(expl_c_res, good_c_res, pop_c_res), Generation, Rule))\n", - "rename!(cval_df, [Symbol(\"$i\") for i in [\"y\", \"ymin\", \"ymax\", \"Generation\", \"Rule\"]]);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(cval_df, x=:Generation, y=:y, ymin=:ymin, ymax=:ymax, color=:Rule, Geom.line, Geom.ribbon,\n", - " Guide.ylabel(\"Average bonus value\"),\n", - " Coord.cartesian(xmin=1, xmax=101),\n", - " Scale.color_discrete_manual(gen_brew_colors(4)[2:4]...),\n", - " style(line_width=2pt, lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2),\n", - " minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now the same for the fitness of the agents." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "full_results_f = Array{Float32, 3}(undef, 200, 101, 100)\n", - "\n", - "for i in 1:100\n", - " full_results_f[:, :, i] = readdlm(\"data/weib_fit$i.txt\")\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bayes_f_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_f[:, k, i][(full_results_type[:, k, i] .== 1)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " bayes_f_res[k, :] = length(unique(cr)) > 1 ? vcat(mean(cr), [confint(OneSampleTTest(cr))...]) : [NaN, NaN, NaN]\n", - "end\n", - "\n", - "expl_f_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_f[:, k, i][(full_results_type[:, k, i] .== 2)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " expl_f_res[k, :] = length(unique(cr)) > 1 ? vcat(mean(cr), [confint(OneSampleTTest(cr))...]) : [NaN, NaN, NaN]\n", - "end\n", - "\n", - "good_f_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_f[:, k, i][(full_results_type[:, k, i] .== 3)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " good_f_res[k, :] = length(unique(cr)) > 1 ? vcat(mean(cr), [confint(OneSampleTTest(cr))...]) : [NaN, NaN, NaN]\n", - "end\n", - "\n", - "pop_f_res = Array{Float64, 2}(undef, 101, 3)\n", - "\n", - "for k in 1:101\n", - " cr = Float64[]\n", - " for j in 1:100\n", - " push!(cr, mean([ full_results_f[:, k, i][(full_results_type[:, k, i] .== 4)] for i in 1:100 ][j]))\n", - " end\n", - " cr = cr[.!isnan.(cr)]\n", - " pop_f_res[k, :] = length(unique(cr)) > 1 ? vcat(mean(cr), [confint(OneSampleTTest(cr))...]) : [NaN, NaN, NaN]\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(\n", - " layer(x=1:101,\n", - " y=expl_f_res[1:101, 1],\n", - " Geom.line,\n", - " Theme(default_color=colorant\"#fc8d62\", line_width=2pt)\n", - " ),\n", - " layer(x=1:101,\n", - " y=good_f_res[1:101, 1],\n", - " Geom.line,\n", - " Theme(default_color=colorant\"#8da0cb\", line_width=2pt)\n", - " ),\n", - " layer(x=1:101,\n", - " y=pop_f_res[1:101, 1],\n", - " Geom.line,\n", - " Theme(default_color=colorant\"#e78ac3\", line_width=2pt)\n", - " ),\n", - " layer(x=1:length(bayes_f_res[.!isnan.(bayes_f_res[:, 1]), 1]),\n", - " y=bayes_f_res[.!isnan.(bayes_f_res[:, 1]), 1],\n", - " ymin=bayes_f_res[.!isnan.(bayes_f_res[:, 1]), 2],\n", - " ymax=bayes_f_res[.!isnan.(bayes_f_res[:, 1]), 3],\n", - " Geom.line,\n", - " Geom.ribbon,\n", - " Theme(default_color=colorant\"#66c2a5\", line_width=2pt, lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2))\n", - " ),\n", - " layer(x=1:101,\n", - " y=expl_f_res[1:101, 1],\n", - " ymin=expl_f_res[1:101, 2],\n", - " ymax=expl_f_res[1:101, 3],\n", - " Geom.line,\n", - " Geom.ribbon,\n", - " Theme(default_color=colorant\"#fc8d62\", lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2))\n", - " ),\n", - " layer(x=1:101,\n", - " y=good_f_res[1:101, 1],\n", - " ymin=good_f_res[1:101, 2],\n", - " ymax=good_f_res[1:101, 3],\n", - " Geom.line,\n", - " Geom.ribbon,\n", - " Theme(default_color=colorant\"#8da0cb\", lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2))\n", - " ),\n", - " layer(x=1:101,\n", - " y=pop_f_res[1:101, 1],\n", - " ymin=pop_f_res[1:101, 2],\n", - " ymax=pop_f_res[1:101, 3],\n", - " Geom.line,\n", - " Geom.ribbon,\n", - " Theme(default_color=colorant\"#e78ac3\", lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2))\n", - " ),\n", - " Guide.xlabel(\"Generation\"),\n", - " Guide.ylabel(\"Average fitness\"),\n", - " Coord.cartesian(xmax=101, ymin=.81),\n", - " Guide.manual_color_key(\"Rule\", [\"Bayes\", \"EXPL\", \"Good\", \"Popper\"], [\"#66c2a5\", \"#fc8d62\", \"#8da0cb\", \"#e78ac3\"]),\n", - " style(minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Summary information about last generations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "last_gen_types = sum([ sum((full_results_type[:, 101, i] .== 1)) for i in 1:100 ]), sum([ sum((full_results_type[:, 101, i] .== 2)) for i in 1:100 ]), sum([ sum((full_results_type[:, 101, i] .== 3)) for i in 1:100 ]), sum([ sum((full_results_type[:, 101, i] .== 4)) for i in 1:100 ])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(x = [\"Bayes\", \"EXPL\", \"Good\", \"Popper\"], y = [last_gen_types...], Geom.bar,\n", - " Guide.xlabel(\"Rule\"),\n", - " Guide.ylabel(\"Count\"),\n", - " Scale.x_discrete,\n", - " style(default_color=colorant\"#66c2a5\", minor_label_font_size=11pt, major_label_font_size=15pt,\n", - " bar_spacing=35pt))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lastE = Float64[]\n", - "for j in 1:100\n", - " append!(lastE, [ full_results_c[:, 101, i][(full_results_type[:, 101, i] .== 2)] for i in 1:100 ][j])\n", - "end\n", - "\n", - "lastG = Float64[]\n", - "for j in 1:100\n", - " append!(lastG, [ full_results_c[:, 101, i][(full_results_type[:, 101, i] .== 3)] for i in 1:100 ][j])\n", - "end\n", - "\n", - "lastP = Float64[]\n", - "for j in 1:100\n", - " append!(lastP, [ full_results_c[:, 101, i][(full_results_type[:, 101, i] .== 4)] for i in 1:100 ][j])\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df = DataFrame(Rule = vcat(fill(\"EXPL\", length(lastE)), fill(\"Good\", length(lastG)), fill(\"Popper\", length(lastP))), C_value = vcat(lastE, lastG, lastP));" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(df, x=:C_value, color=:Rule, Geom.density(bandwidth=.005),\n", - " Coord.cartesian(xmin=-.0075, xmax=.261),\n", - " Scale.color_discrete_manual(gen_brew_colors(4)[2:4]...),\n", - " Guide.xlabel(\"Bonus value\"),\n", - " Guide.ylabel(\"Density\"),\n", - " style(line_width=2.65pt, minor_label_font_size=10pt, major_label_font_size=14pt,\n", - " key_label_font_size=11pt, key_title_font_size=13pt,\n", - " colorkey_swatch_shape=:square))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Statistics" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using RCall" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lastEF = Float64[]\n", - "for j in 1:100\n", - " append!(lastEF, [ full_results_f[:, 101, i][(full_results_type[:, 101, i] .== 2)] for i in 1:100 ][j])\n", - "end\n", - "\n", - "lastGF = Float64[]\n", - "for j in 1:100\n", - " append!(lastGF, [ full_results_f[:, 101, i][(full_results_type[:, 101, i] .== 3)] for i in 1:100 ][j])\n", - "end\n", - "\n", - "lastPF = Float64[]\n", - "for j in 1:100\n", - " append!(lastPF, [ full_results_f[:, 101, i][(full_results_type[:, 101, i] .== 4)] for i in 1:100 ][j])\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x = vcat(lastEF, lastGF, lastPF)\n", - "y = vcat(fill(\"EXPL\", length(lastEF)), fill(\"Good\", length(lastGF)), fill(\"Popper\", length(lastPF)))\n", - "R\"summary(aov($x ~ $y))\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "R\"pairwise.t.test($x, $y, p.adj = 'bonf')\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "R\"library(lsr)\"\n", - "R\"eSq <- etaSquared(aov($x ~ $y))[1]\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "R\"library(emmeans)\"\n", - "R\"df <- data.frame(x = $x, y = $y)\"\n", - "R\"mod <- lm(x ~ y, data = df)\"\n", - "R\"emmeans(mod, ~ y)\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x1 = vcat(lastE, lastG, lastP)\n", - "y1 = vcat(fill(\"EXPL\", length(lastE)), fill(\"Good\", length(lastG)), fill(\"Popper\", length(lastP)))\n", - "R\"df <- data.frame(x = $x1, y = $y1)\"\n", - "R\"mod <- lm(x ~ y, data = df)\"\n", - "R\"emmeans(mod, ~ y)\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.3.1", - "language": "julia", - "name": "julia-1.3" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.3.1" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -}