{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# S1 - Sample-based implementation of Blahut-Arimoto iteration\n", "This notebook is part of the supplementary material for: \n", "Genewein T., Leibfried F., Grau-Moya J., Braun D.A. (2015) *Bounded rationality, abstraction and hierarchical decision-making: an information-theoretic optimality principle*, Frontiers in Robotics and AI. \n", "\n", "More information on how to run the notebook on the accompanying [github repsitory](https://github.com/tgenewein/BoundedRationalityAbstractionAndHierarchicalDecisionMaking) where you can also find updated versions of the code and notebooks.\n", "\n", "This notebook in mentioned in Section 2.3. Due to time- and space-limitations the results of this notebook are not in the paper." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Disclaimer\n", "This notebook provides a proof-of-concept implementation of a naive sample-based Blahut-Arimoto iteration scheme. Neither the code nor the notebook have been particularly polished or tested. There is a short theory-bit in the beginning of the notebook but most of the explanations are brief and mixed into the code as comments." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Free energy rejection sampling\n", "The solution to a free energy variational problem (Section 2.1 in the paper) has the form of a Boltzmann distribution\n", "$$p(y) = \\frac{1}{Z}p_0(y)e^{\\beta U(y)},$$\n", "where $Z=\\sum_y p_0(y)e^{\\beta U(y)}$ denotes the partition sum, $p_0(y)$ is a prior distribution and $U(y)$ is the utility function. The inverse temperature $\\beta$ can be interpreted as a resource parameter and it governs how far the posterior $p(y)$ can deviate from the prior (measured as a KL-divergence) - see the paper Section 2.1 for details.\n", "\n", "For a decision-maker, it suffices to obtain a sample from $p(y)$ and act according to that sample, rather than computing the full distribution $p(y)$. A simple scheme to sample from $p(y)$ is given by rejection sampling.\n", "\n", "**Rejection sampling** \n", "Goal: get a sample $y$ from the distribution $f(y)$. Draw from a uniform distribution $u\\sim \\mathcal{U}(0,1)$ and from a proposal distribution $y\\sim g(y)$. If $u < \\frac{f(y)}{M g(y)}$, accept the sample as a sample from $f(y)$, otherwise reject the sample and repeat. $M$ is a constant that ensures that $M g(y) \\geq f(y)~\\forall y$. Note that rejection sampling also works for sampling from an unnormalized distribution as long as $M$ is chosen accordingly.\n", "\n", "For the free-energy problem, we want a sample from $p(y)\\propto f(y) = p_0(y)e^{\\beta U(y)}$. We choose $g(y)=p_0(y)$ and set $M=e^{\\beta U_{max}}$, where $U_{max}=\\underset{y}{max}~U(y)$\n", "\n", "**Finally we get the following rejection sampling scheme:** \n", "* draw from a uniform distribution $u\\sim \\mathcal{U}(0,1)$ \n", "* draw from the proposal distribution $x\\sim p_0(y)$ (*the prior*)\n", " * if $u < \\frac{\\exp(\\beta U(y))}{\\exp(\\beta U_\\mathrm{max})}$ accept the sample as a sample from the posterior $p(y)$\n", " * otherwise reject the sample (and re-sample). \n", "\n", "\n", "## Rate distortion rejection sampling\n", "The solution to the rate distortion problem looks very similar to the Boltzmann distribution in the free-energy case. However, there is one crucial difference: in the free-energy case, the prior is an arbitrary distribution - in the rate distortion case, the prior is replaced by the marginal distribution, which leads to a set of self-consistent equations\n", "$$\\begin{align}\n", "p^*(a|w)&=\\frac{1}{Z(w)}p(a)e^{\\beta U(a,w)} \\\\\n", "p(a)&=\\sum_w p(w)p(a|w)\n", "\\end{align}$$\n", "*After* convergence of the Blahut-Arimoto iterations, the marginal $p(a)$ can just be treated like a prior and the rejection sampling scheme described above can straightforwardly be used. However, when initializing with an arbitary marginal distribution $\\hat{p}(a)$ the iterations must be performed in a sample-based manner until convergence.\n", "\n", "Here, we do this in a naive and very straightforward way: we represent $\\hat{p}(a)$ simply through counters (a categorical distribution). Then we do the following:\n", "1. Draw a number of samples (a batch) from $\\hat{p}^*(a|w)=\\frac{1}{Z(w)}\\hat{p}(a)e^{\\beta U(a,w)}$ using the rejection sampling scheme.\n", "2. Update $\\hat{p}(a)$ with the accepted samples obtained in step 1. There are different possibilities for the update step\n", " 1. Simply increase the counters for each accepted $a$ and re-normalize (no-forgetting)\n", " 2. Reset the counters for $a$ and use only the last batch of accepted samples to empirically estimate $p(a)$ (full-forgetting)\n", " 3. Use an exponentially decaying window over the last samples to update the empirical estimate of $p(a)$ (not implemented in this notebook).\n", " 4. Use a parametric model for $p_theta(a)$ and then perform some moment-matching or use a gradient-based update rule to adjust the parameters $\\theta$ (not implemented in this notebook).\n", "3. Repeat until convergence (or here for simplicity: for a fixed number of steps)\n", "\n", "Additionally, this notebook allows for some burn-in time, where after a certain number of iterations of 1. and 2. (i.e. after the \"burn-in\") the counters for $\\hat{p}(a)$ are reset. This naive scheme seems to work but it is unclear how to choose the batch-size (number of samples from $\\hat{p}^*(a|w)$ to take before performing an update step on $\\hat{p}(a)$), how to set the burn-in phase, etc.\n", "\n", "\n", "In the notebook below, you can try different batch-sizes and different burn-in times and you can compare full-forgetting against no forgetting (i.e. no resetting of the counters). " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#only run this once\n", "include(\"RateDistortionDecisionMaking.jl\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load the taxonomy example as a testbed" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#load taxonomy example\n", "using RateDistortionDecisionMaking, DataFrames, Gadfly, Distributions\n", "\n", "#set up taxonomy example\n", "include(\"TaxonomyExample.jl\")\n", "w_vec, w_strings, a_vec, a_strings, p_w, U = setuptaxonomy()\n", "\n", "#pre-compute utilities, find maxima\n", "U_pre, Umax = setuputilityarrays(a_vec,w_vec,U)\n", "\n", "#initialize p(a) uniformly\n", "num_acts = length(a_vec)\n", "pa_init = ones(num_acts)/num_acts;\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Set up functions for sampling and run on the example from above\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#Performs rejection sampling with a constant (scaled uniform) envelope\n", "#using a softmax acceptance-rejection criterion.\n", "#prop_dist .. proposal distribution (must be an instance of type ::Distribution)\n", "#nsamps ..... desired number of samples (scalar)\n", "#maxsteps ... maximum number of acceptance-rejection steps (scalar, must be ≧ nsamps)\n", "#β .......... softmax parameter\n", "#lh ......... likelihood value (vector of length N)\n", "#maxlh ...... maximum value that the likelihood can take (scalar)\n", "function rej_samp_const(prop_dist::Distribution, nsamps::Integer, maxsteps::Integer, β::Number, lh::Vector, maxlh::Number) \n", " #initialize\n", " samps = zeros(nsamps)\n", " acc_cnt = 0 #acceptance-counter\n", " if(maxsteps < nsamps)\n", " maxsteps = nsamps\n", " end\n", " \n", " k=0 #use this to make sure that k is still available after the loop\n", " for k in 1:maxsteps\n", " u=rand(1) #sample from uniform between (0,1)\n", " index = rand(prop_dist) #sample from proposal\n", " \n", " ratio = exp(β*lh[index])/exp(β*maxlh)\n", " if u[1]= can not handle arrays\n", " #if we enter here, accept the sample \n", " acc_cnt = acc_cnt + 1 \n", " samps[acc_cnt] = index\n", " \n", " if(acc_cnt == nsamps)\n", " #we have enough samples, exit loop\n", " break\n", " end\n", " end\n", " end\n", " \n", " if(k==maxsteps)\n", " warn(\"[RejSampConst] Maximum number of steps reached - number of samples is potentially lower than nsamps!\\n\")\n", " end\n", " \n", " #store all accepted samples (this can be less than nsamps if maxsteps is too low or acceptance-rate is low)\n", " samples = samps[1:acc_cnt]\n", " \n", " #compute acceptance ratio\n", " acc_ratio = acc_cnt/k\n", " \n", " return samples, acc_ratio\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#marginal is simply represented by counters (i.e. by frequencies)\n", "function init_marginal_representation_ctrs(pa_init::Vector)\n", " return pa_init\n", "end\n", "\n", "#this updates the marginal over actions p(a) using a counter-representation\n", "#this function counts the number of times each action-index occurs in sampled_indices\n", "#these counts are then added to the current marginal_ctrs. Optionally, the counters are reset\n", "#before adding the new samples (=hard forgetting).\n", "function update_marginal_ctrs(sampled_indices::Vector, marginal_ctrs::Vector; reset_ctrs::Bool=false) \n", " #TODO: perhaps replace hard-resetting with an exponential decay?\n", " \n", " p_ctrs = marginal_ctrs\n", " card_p = length(p_ctrs)\n", " \n", " #reset counters for marginal? (make sure every entry is non-zero!)\n", " if reset_ctrs\n", " p_ctrs = ones(card_p)/card_p\n", " end\n", "\n", " #update marginal counters using a histogram to do the counting (bin-borders have to be set manually!)\n", " e,p_counts = hist(sampled_indices,0.5:1:(card_p+0.5)) \n", " p_ctrs = p_ctrs + p_counts\n", " \n", " #normalize to get the updated marginal\n", " p_sampled = p_ctrs / sum(p_ctrs)\n", " \n", " return p_sampled, p_ctrs #return the probability-vector, but also the representation of the marginal (as counts)\n", "end\n", "\n", "\n", "\n", "#function for BA sampling\n", "#burnin_ratio specifies the ratio of outer iterations that will not count\n", "#towards computation of the final marginal distribution (counters will be blocked)\n", "#reset_marginal_ctrs specifies whether the marginal is computed with the samples of the last\n", "#iteration only (=hard forgetting by resetting counters) or whether the marginal\n", "#is computed with all samples of all iterations (=no forgetting)\n", "function BAsampling(pa_init::Vector, β::Number, U_pre::Matrix, Umax::Vector, pw::Vector, \n", " nsteps_marginalupdate::Integer, nsteps_conditionalupdate::Integer;\n", " burnin_ratio::Real=0.7, max_rejsamp_steps::Integer=200,\n", " compute_performance::Bool=false, performance_as_dataframe::Bool=false,\n", " performance_per_iteration::Bool=false,\n", " init_marg_func::Function=init_marginal_representation_ctrs,\n", " update_marg_func::Function=update_marginal_ctrs, update_func_args...)\n", " \n", " #compute cardinality, check size of U_pre\n", " card_a = length(pa_init)\n", " card_w = length(pw)\n", " if size(U_pre) != (card_a, card_w)\n", " error(\"Size mismatch of U_pre and pa_init or pw!\")\n", " end\n", " \n", " #check that burnin_ratio is really a ratio\n", " if (burnin_ratio < 0) || (burnin_ratio > 1)\n", " error(\"burnin_ratio must be a number between 0 and 1.\")\n", " end\n", " \n", " #if performance measures don't need to be returned, don't compute them per iteration\n", " if compute_performance==false\n", " performance_per_iteration = false\n", " end \n", " #preallocate if necessary\n", " if performance_per_iteration \n", " I_i = zeros(maxiter)\n", " Ha_i = zeros(maxiter)\n", " Hagw_i = zeros(maxiter)\n", " EU_i = zeros(maxiter)\n", " RDobj_i = zeros(maxiter)\n", " end\n", " \n", " #initialize sampling distributions\n", " pw_dist = Categorical(pw) #proposal distribution \n", " pagw_ctrs = ones(card_a, card_w) #counters for conditional distribution \n", " pa_sampled = pa_init #marginal distribution\n", " \n", " #initialize the marginal representation\n", " pa_ctrs = init_marg_func(pa_init)\n", "\n", " burnin_triggered=false\n", " #outer loop - in each iteration the marginal is updated\n", " iter=0\n", " for iter in 1:nsteps_marginalupdate \n", " a_samples = zeros(nsteps_conditionalupdate) #this will hold the samples from p(a|w) during inner loop\n", " \n", " #inner loop - in each step a sample is drawn from the conditional and stored for\n", " #for the batch-update of the marginal\n", " for j in 1:nsteps_conditionalupdate\n", " #draw a w sample\n", " w_samp = rand(pw_dist)\n", "\n", " #draw a sample from p(a|w) using the current estimate of p(a) as proposal distribution using rejection sampling\n", " agw_samp, acc_ratio = rej_samp_const(Categorical(pa_sampled), 1, max_rejsamp_steps, β, U_pre[:,w_samp], Umax[w_samp])\n", " a_samples[j] = agw_samp[1]\n", " \n", "\n", " #update conditional counters\n", " pagw_ctrs[agw_samp, w_samp] += 1\n", " end\n", " \n", " #very simple burn-in: simply reset counters\n", " if (iter >(nsteps_marginalupdate)*burnin_ratio) && (!burnin_triggered)\n", " burnin_triggered = true\n", " pagw_ctrs = ones(card_a, card_w) \n", " end\n", "\n", " #update marginal with samples drawn in inner loop\n", " pa_sampled, pa_ctrs = update_marg_func(a_samples, pa_ctrs; update_func_args...) \n", " \n", " \n", " #compute entropic quantities (if requested with additional parameter)\n", " if performance_per_iteration\n", " #compute sample-based conditional p(a|w)\n", " pagw_sampled = zeros(card_a, card_w)\n", " for i in 1:card_w\n", " pagw_sampled[:,i] = pagw_ctrs[:,i] / sum(pagw_ctrs[:,i])\n", " end\n", " I_i[iter], Ha_i[iter], Hagw_i[iter], EU_i[iter], RDobj_i[iter] = analyzeBAsolution(pw, pa_sampled, pagw_sampled, U_pre, β)\n", " end\n", " end\n", "\n", " #compute conditionals using the sample-counts of the previous inner loops\n", " #the burn-in parameter specifies how many of the inner loops are discarded\n", " pagw_sampled = zeros(card_a, card_w)\n", " for i in 1:card_w\n", " pagw_sampled[:,i] = pagw_ctrs[:,i] / sum(pagw_ctrs[:,i])\n", " end\n", "\n", " #return results\n", " if compute_performance == false\n", " return pagw_sampled, pa_sampled\n", " else \n", " if performance_per_iteration == false\n", " #compute performance measures for final solution\n", " I, Ha, Hagw, EU, RDobj = analyzeBAsolution(pw, pa_sampled, pagw_sampled, U_pre, β)\n", " else\n", " #\"cut\" valid results from preallocated vector\n", " I = I_i[1:iter]\n", " Ha = Ha_i[1:iter]\n", " Hagw = Hagw_i[1:iter]\n", " EU = EU_i[1:iter]\n", " RDobj = RDobj_i[1:iter]\n", " end\n", "\n", " #if needed, transform to data frame\n", " if performance_as_dataframe == false\n", " return pagw_sampled, pa_sampled, I, Ha, Hagw, EU, RDobj\n", " else\n", " performance_df = performancemeasures2DataFrame(I, Ha, Hagw, EU, RDobj)\n", " return pagw_sampled, pa_sampled, performance_df \n", " end\n", " end\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#example call and also plot evolution of performance measueres\n", "maxiter = 10000\n", "β = 1.2\n", "nsteps_marg = 500\n", "nsteps_cond = 750\n", "pagw_s,pa_s,perf = BAsampling(pa_init, β, U_pre, Umax, p_w, nsteps_marg, nsteps_cond,\n", " burnin_ratio=0.7, max_rejsamp_steps=500, reset_ctrs=false,\n", " compute_performance=true, performance_as_dataframe=true, performance_per_iteration=true)\n", "\n", "plt_cond = visualizeBAconditional(pagw_s,a_vec,w_vec,a_strings,w_strings)\n", "\n", "#instead of using a range of β-values (as for the standard-performance plot), \n", "#use a vector indicating the iteration\n", "niter = size(perf,1)\n", "plt_perf_entropy, plt_perf_utility, plt_rateutility = plotperformancemeasures(perf,collect(1:niter),\n", " suppress_vis=true, xlabel_perf=\"Iteration\")\n", "\n", "#TODO: somehow the \"Iteration\" label above doesn't seem to work!\n", "\n", "display(vstack(plt_perf_entropy, plt_perf_utility))\n", "display(plt_cond)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## [Interact] Change the parameters in the code-cell above to explore the sampling scheme and its solutions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare sampling-solutions against analytical results\n", "\n", "Below, we will average over several sampling-runs at different temperatures $\\beta$ to see a difference between the analytical solutions and the sample-based solutions (with and without forgetting)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#compute theoretical result for rate-dutility curve\n", "ε = 0.0001 #convergence critetion for BAiterations\n", "maxiter = 10000\n", "β_sweep = collect(0.01:0.05:3)\n", "nβ = length(β_sweep)\n", "\n", "#preallocate\n", "I = zeros(nβ)\n", "Ha = zeros(nβ)\n", "Hagw = zeros(nβ)\n", "EU = zeros(nβ)\n", "RDobj = zeros(nβ)\n", "\n", "#sweep through β values and perfomr Blahut-Arimoto iterations for each value\n", "for i=1:nβ \n", " pagw, pa, I[i], Ha[i], Hagw[i], EU[i], RDobj[i] = BAiterations(pa_init, β_sweep[i], U_pre, p_w, ε, maxiter,compute_performance=true) \n", "end\n", "\n", "#show rate-utility curve (shaded region is theoretically infeasible)\n", "perf_res_analytical = performancemeasures2DataFrame(I, Ha, Hagw, EU, RDobj); \n", "plot_perf_entropy, plot_perf_util, plot_rateutility = plotperformancemeasures(perf_res_analytical, β_sweep, suppress_vis=true);\n", "display(plot_rateutility)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#run the smapling for different temperatures and repeat each run n-times\n", "#then plot these results against the rate-utility curve (based on the closed-form solutions)\n", "βrange_samp = [0.1, 0.25, 0.5, 0.8, 1.2, 1.4, 1.6, 2]\n", "#βrange_samp = [1.2, 2]\n", "\n", "nruns = 10; #number of runs per β point\n", "\n", "nsteps_marg = 500\n", "nsteps_cond = 750\n", "burnin_ratio = 0.8\n", "max_rejsamp_steps=500 #maximum number of steps for sampling from the conditional\n", "\n", "\n", "nconditions = size(βrange_samp,1)*nruns\n", "I_sampled = zeros(2*nconditions)\n", "Ha_sampled = zeros(2*nconditions)\n", "EU_sampled = zeros(2*nconditions)\n", "βval = zeros(2*nconditions)\n", "ResetCtrs = falses(2*nconditions)\n", "\n", "#first run with reset_ctrs to false\n", "reset_ctrs = false #if true, the ctrs for the marginal are reset in each iteration (=\"hard\" forgetting)\n", "for b in 1:nconditions\n", " println(\"BA Sampling, run $b of $(2*nconditions)\")\n", " β = βrange_samp[round(Int,ceil(b/nruns))]\n", " \n", " pagw_s, pa_s, I, Ha, Hagw, EU, RDobj = BAsampling(pa_init, β, U_pre, Umax, p_w, nsteps_marg, nsteps_cond,\n", " reset_ctrs=reset_ctrs, burnin_ratio=burnin_ratio,\n", " max_rejsamp_steps=max_rejsamp_steps, compute_performance=true)\n", " \n", " I_sampled[b] = I\n", " Ha_sampled[b] = Ha\n", " EU_sampled[b] = EU\n", " βval[b] = β\n", " ResetCtrs[b] = reset_ctrs\n", "end\n", "\n", "#second run with reset_ctrs to true\n", "reset_ctrs = true #if true, the ctrs for the marginal are reset in each iteration (=\"hard\" forgetting)\n", "for b in 1:nconditions\n", " println(\"BA Sampling, run $(nconditions+b) of $(2*nconditions)\")\n", " β = βrange_samp[round(Int,ceil(b/nruns))]\n", " \n", " pagw_s, pa_s, I, Ha, Hagw, EU, RDobj = BAsampling(pa_init, β, U_pre, Umax, p_w, nsteps_marg, nsteps_cond,\n", " reset_ctrs=reset_ctrs, burnin_ratio=burnin_ratio,\n", " max_rejsamp_steps=max_rejsamp_steps, compute_performance=true)\n", " \n", " I_sampled[nconditions+b] = I\n", " Ha_sampled[nconditions+b] = Ha\n", " EU_sampled[nconditions+b] = EU\n", " βval[nconditions+b] = β\n", " ResetCtrs[nconditions+b] = reset_ctrs\n", "end\n", "\n", "#wrap data in DataFrame for convenient plotting\n", "res_sampled = DataFrame(β=βval, I_aw=I_sampled, H_a=Ha_sampled, E_U=EU_sampled, Forgetting=ResetCtrs)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#compute theoretical result for same set of temperatures\n", "ε = 0.0001 #convergence critetion for BAiterations\n", "maxiter = 10000\n", "nβ = length(βrange_samp)\n", "\n", "#preallocate\n", "I = zeros(nβ)\n", "Ha = zeros(nβ)\n", "Hagw = zeros(nβ)\n", "EU = zeros(nβ)\n", "RDobj = zeros(nβ)\n", "\n", "#sweep through β values and perfomr Blahut-Arimoto iterations for each value\n", "for i=1:nβ \n", " pagw, pa, I[i], Ha[i], Hagw[i], EU[i], RDobj[i] = BAiterations(pa_init, βrange_samp[i], U_pre, p_w, ε, maxiter,compute_performance=true) \n", "end\n", "\n", "#show rate-utility curve (shaded region is theoretically infeasible)\n", "perf_res_analytical_samp = performancemeasures2DataFrame(I, Ha, Hagw, EU, RDobj)\n", "perf_res_analytical_samp[:β] = βrange_samp;" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#plot the solutions from the sampling runs into the (analytical) rate-utility plot\n", "plot(Guide.ylabel(\"E[U]\"), Guide.xlabel(\"I(A;W) [bits]\"),\n", " Guide.title(\"Rate-Utility curve (dots show sampling solutions)\"), BAtheme(), BAdiscretecolorscale(2),\n", " layer(res_sampled,y=\"E_U\",x=\"I_aw\",Geom.point,color=\"Forgetting\"),\n", " layer(perf_res_analytical_samp,y=\"E_U\",x=\"I_aw\",Geom.point),\n", " layer(perf_res_analytical,y=\"E_U\",x=\"I_aw\",Geom.line),\n", " layer(perf_res_analytical,y=\"E_U\",x=\"I_aw\",ymin=\"E_U\",ymax=ones(size(perf_res_analytical,1))*maximum(perf_res_analytical[:E_U]),\n", " Geom.ribbon,BAtheme(default_color=colorant\"green\"))\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#compute the mean for each group of points (that were produced with the same beta and same forgetting setting)\n", "res_samp_aggregated = aggregate(res_sampled, [:β,:Forgetting], mean)\n", "\n", "#plot the mean-solutions from the sampling runs into the (analytical) rate-utility plot\n", "plt_samp = plot(Guide.ylabel(\"E[U]\"), Guide.xlabel(\"I(A;W) [bits]\"),\n", " Guide.title(\"Rate-Utility curve (dots show mean sampling solutions)\"), BAtheme(), BAdiscretecolorscale(2),\n", " layer(res_samp_aggregated,y=\"E_U_mean\",x=\"I_aw_mean\",Geom.point,Geom.line(preserve_order=true),color=\"Forgetting\"),\n", " layer(perf_res_analytical_samp,y=\"E_U\",x=\"I_aw\",Geom.point),\n", " layer(perf_res_analytical,y=\"E_U\",x=\"I_aw\",Geom.line),\n", " layer(perf_res_analytical,y=\"E_U\",x=\"I_aw\",ymin=\"E_U\",ymax=ones(size(perf_res_analytical,1))*maximum(perf_res_analytical[:E_U]),\n", " Geom.ribbon,BAtheme(default_color=colorant\"green\"))\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#store the plots\n", "#draw(SVG(\"Figures/RateUtilityCurve.svg\", 8.5cm, 7cm), plot_rateutility)\n", "\n", "#draw(SVG(\"Figures/SamplingCond.svg\", 8.5cm, 9cm), plt_cond)\n", "#draw(SVG(\"Figures/BASampling.svg\", 13cm, 11cm), plt_samp)\n", "\n", "#plot_samp = vstack(plt_samp, plt_cond)\n", "#draw(SVG(\"Figures/SamplingResults.svg\", 18cm,16cm),plot_samp)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#try changing the number of inner-/outer-loop iterations\n", "#try changing the burn-in ratio\n", "#try soft-forgetting (with exponential-decay window)\n", "\n", "#forgetting seems to do better than no forgetting (in terms of being closer to the rate-utility curve),\n", "#but it also seems that the points with forgetting tend to have a lower I(A;O) (and also a lower E[U]) - even though\n", "#the temperatures are the same." ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.4.0", "language": "julia", "name": "julia-0.4" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.4.0" } }, "nbformat": 4, "nbformat_minor": 0 }