docs/source/tuto_network_calibration/clustering_dhist.ipynb

   1 {
   2  "cells": [
   3   {
   4    "cell_type": "markdown",
   5    "metadata": {},
   6    "source": [
   7     "# Clustering data using Dhist"
   8    ]
   9   },
  10   {
  11    "cell_type": "markdown",
  12    "metadata": {},
  13    "source": [
  14     "Use diagonally cut histogram to model the variability.\n",
  15     "\n",
  16     "Dhist creates a non-uniform histogram which is sampled later to estimate the intercept."
  17    ]
  18   },
  19   {
  20    "cell_type": "code",
  21    "execution_count": 3,
  22    "metadata": {
  23     "tags": [
  24      "parameters"
  25     ]
  26    },
  27    "outputs": [],
  28    "source": [
  29     "seed=42\n",
  30     "df_str=\"\"\n",
  31     "metric=\"\"\n",
  32     "coeff_col_name=\"\"\n",
  33     "coeff_value=0\n",
  34     "log=TRUE\n",
  35     "xbr_file=\"/tmp/dhist_xbr.csv\"\n",
  36     "height_file=\"/tmp/dhist_height.csv\""
  37    ]
  38   },
  39   {
  40    "cell_type": "markdown",
  41    "metadata": {},
  42    "source": [
  43     "We use a special method to create non-uniform histograms to represent the noise in IO operations.\n",
  44     "\n",
  45     "Unable to install the library properly, I copied the important methods here.\n",
  46     "\n",
  47     "Copied from: https://rdrr.io/github/dlebauer/pecan-priors/src/R/plots.R"
  48    ]
  49   },
  50   {
  51    "cell_type": "code",
  52    "execution_count": 5,
  53    "metadata": {},
  54    "outputs": [],
  55    "source": [
  56     "set.seed(seed) #forcing seed for reproductibility\n",
  57     "\n",
  58     "#' Variable-width (dagonally cut) histogram\n",
  59     "#'\n",
  60     "#'\n",
  61     "#' When constructing a histogram, it is common to make all bars the same width.\n",
  62     "#' One could also choose to make them all have the same area.\n",
  63     "#' These two options have complementary strengths and weaknesses; the equal-width histogram oversmooths in regions of high density, and is poor at identifying sharp peaks; the equal-area histogram oversmooths in regions of low density, and so does not identify outliers.\n",
  64     "#' We describe a compromise approach which avoids both of these defects. We regard the histogram as an exploratory device, rather than as an estimate of a density.\n",
  65     "#' @title Diagonally Cut Histogram\n",
  66     "#' @param x is a numeric vector (the data)\n",
  67     "#' @param a is the scaling factor, default is 5 * IQR\n",
  68     "#' @param nbins is the number of bins, default is assigned by the Stuges method\n",
  69     "#' @param rx  is the range used for the left of the left-most bin to the right of the right-most bin\n",
  70     "#' @param eps used to set artificial bound on min width / max height of bins as described in Denby and Mallows (2009) on page 24.\n",
  71     "#' @param xlab is label for the x axis\n",
  72     "#' @param plot = TRUE produces the plot, FALSE returns the heights, breaks and counts\n",
  73     "#' @param lab.spikes = TRUE labels the \\% of data in the spikes\n",
  74     "#' @return list with two elements, heights of length n and breaks of length n+1 indicating the heights and break points of the histogram bars.\n",
  75     "#' @author Lorraine Denby, Colin Mallows\n",
  76     "#' @references Lorraine Denby, Colin Mallows. Journal of Computational and Graphical Statistics. March 1, 2009, 18(1): 21-31. doi:10.1198/jcgs.2009.0002.\n",
  77     "dhist<-function(x, a=5*iqr(x),\n",
  78     "                 nbins=nclass.Sturges(x), rx = range(x,na.rm = TRUE),\n",
  79     "                 eps=.15, xlab = \"x\", plot = TRUE,lab.spikes = TRUE)\n",
  80     " {\n",
  81     "\n",
  82     "   if(is.character(nbins))\n",
  83     "     nbins <- switch(casefold(nbins),\n",
  84     "                     sturges = nclass.Sturges(x),\n",
  85     "                     fd = nclass.FD(x),\n",
  86     "                     scott = nclass.scott(x),\n",
  87     "                     stop(\"Nclass method not recognized\"))\n",
  88     "   else if(is.function(nbins))\n",
  89     "     nbins <- nbins(x)\n",
  90     "\n",
  91     "   x <- sort(x[!is.na(x)])\n",
  92     "   if(a == 0)\n",
  93     "     a <- diff(range(x))/100000000\n",
  94     "   if(a != 0 & a != Inf) {\n",
  95     "     n <- length(x)\n",
  96     "     h <- (rx[2] + a - rx[1])/nbins\n",
  97     "     ybr <- rx[1] + h * (0:nbins)\n",
  98     "     yupper <- x + (a * (1:n))/n\n",
  99     "                                         # upper and lower corners in the ecdf\n",
 100     "     ylower <- yupper - a/n\n",
 101     "                                         #\n",
 102     "     cmtx <- cbind(cut(yupper, breaks = ybr), cut(yupper, breaks = \n",
 103     "                                 ybr, left.include = TRUE), cut(ylower, breaks = ybr),\n",
 104     "                   cut(ylower, breaks = ybr, left.include = TRUE))\n",
 105     "     cmtx[1, 3] <- cmtx[1, 4] <- 1\n",
 106     "                                         # to replace NAs when default r is used\n",
 107     "     cmtx[n, 1] <- cmtx[n, 2] <- nbins\n",
 108     "                                         #\n",
 109     "                                         #checksum <- apply(cmtx, 1, sum) %% 4\n",
 110     "     checksum <- (cmtx[, 1] + cmtx[, 2] + cmtx[, 3] + cmtx[, 4]) %%\n",
 111     "     4\n",
 112     "                                         # will be 2 for obs. that straddle two bins\n",
 113     "     straddlers <- (1:n)[checksum == 2]\n",
 114     "                                         # to allow for zero counts\n",
 115     "     if(length(straddlers) > 0) {\n",
 116     "       counts <- table(c(1:nbins, cmtx[ - straddlers, 1]))\n",
 117     "     } else {\n",
 118     "       counts <- table(c(1:nbins, cmtx[, 1]))\n",
 119     "     }\n",
 120     "     counts <- counts - 1\n",
 121     "                                         #\n",
 122     "     if(length(straddlers) > 0) {\n",
 123     "       for(i in straddlers) {\n",
 124     "         binno <- cmtx[i, 1]\n",
 125     "         theta <- ((yupper[i] - ybr[binno]) * n)/a\n",
 126     "         counts[binno - 1] <- counts[binno - 1] + (\n",
 127     "                                                   1 - theta)\n",
 128     "         counts[binno] <- counts[binno] + theta\n",
 129     "       }\n",
 130     "     }\n",
 131     "     xbr <- ybr\n",
 132     "     xbr[-1] <- ybr[-1] - (a * cumsum(counts))/n\n",
 133     "     spike<-eps*diff(rx)/nbins\n",
 134     "     flag.vec<-c(diff(xbr)<spike,F)\n",
 135     "     if ( sum(abs(diff(xbr))<=spike) >1) {\n",
 136     "       xbr.new<-xbr\n",
 137     "       counts.new<-counts\n",
 138     "       diff.xbr<-abs(diff(xbr))\n",
 139     "       amt.spike<-diff.xbr[length(diff.xbr)]\n",
 140     "       for (i in rev(2:length(diff.xbr))) {\n",
 141     "         if (diff.xbr[i-1]<=spike&diff.xbr[i]<=spike&\n",
 142     "             !is.na(diff.xbr[i])) {\n",
 143     "           amt.spike<-amt.spike+diff.xbr[i-1]\n",
 144     "           counts.new[i-1]<-counts.new[i-1]+counts.new[i]\n",
 145     "           xbr.new[i]<-NA\n",
 146     "           counts.new[i]<-NA\n",
 147     "           flag.vec[i-1]<-T\n",
 148     "         }\n",
 149     "         else amt.spike<-diff.xbr[i-1]\n",
 150     "       }\n",
 151     "       flag.vec<-flag.vec[!is.na(xbr.new)]\n",
 152     "       flag.vec<-flag.vec[-length(flag.vec)]\n",
 153     "       counts<-counts.new[!is.na(counts.new)]\n",
 154     "       xbr<-xbr.new[!is.na(xbr.new)]\n",
 155     "\n",
 156     "     }\n",
 157     "     else flag.vec<-flag.vec[-length(flag.vec)]\n",
 158     "     widths <- abs(diff(xbr))\n",
 159     "     ## N.B. argument \"widths\" in barplot must be xbr\n",
 160     "     heights <- counts/widths\n",
 161     "   }\n",
 162     "   bin.size <- length(x)/nbins\n",
 163     "   cut.pt <- unique(c(min(x) - abs(min(x))/1000,\n",
 164     "                      approx(seq(length(x)), x, (1:(nbins - 1)) * bin.size, rule = 2)$y, max(x)))\n",
 165     "   aa <- hist(x, breaks = cut.pt, plot = FALSE, probability = TRUE)\n",
 166     "   if(a == Inf) {\n",
 167     "     heights <- aa$counts\n",
 168     "     xbr <- aa$breaks\n",
 169     "   }\n",
 170     "   amt.height<-3\n",
 171     "   q75<-quantile(heights,.75)\n",
 172     "   if (sum(flag.vec)!=0) {\n",
 173     "     amt<-max(heights[!flag.vec])\n",
 174     "     ylim.height<-amt*amt.height\n",
 175     "     ind.h<-flag.vec&heights> ylim.height\n",
 176     "     flag.vec[heights<ylim.height*(amt.height-1)/amt.height]<-F\n",
 177     "     heights[ind.h] <- ylim.height\n",
 178     "   }\n",
 179     "   amt.txt<-0\n",
 180     "   end.y<-(-10000)\n",
 181     "   if(plot) {\n",
 182     "     barplot(heights, abs(diff(xbr)), space = 0, density = -1, xlab = \n",
 183     "             xlab, plot = TRUE, xaxt = \"n\",yaxt='n')\n",
 184     "     at <- pretty(xbr)\n",
 185     "     axis(1, at = at - xbr[1], labels = as.character(at))\n",
 186     "     if (lab.spikes) {\n",
 187     "       if (sum(flag.vec)>=1) {\n",
 188     "         usr<-par('usr')\n",
 189     "         for ( i in seq(length(xbr)-1)) {\n",
 190     "           if (!flag.vec[i]) {\n",
 191     "             amt.txt<-0\n",
 192     "             if (xbr[i]-xbr[1]<end.y) amt.txt<-1\n",
 193     "           }\n",
 194     "           else {\n",
 195     "             amt.txt<-amt.txt+1\n",
 196     "             end.y<-xbr[i]-xbr[1]+3*par('cxy')[1]\n",
 197     "           }\n",
 198     "           if (flag.vec[i]) {\n",
 199     "             txt<-paste(' ',format(round(counts[i]/\n",
 200     "                                         sum(counts)*100)),'%',sep='')\n",
 201     "             par(xpd = TRUE)\n",
 202     "             text(xbr[i+1]-xbr[1],ylim.height-par('cxy')[2]*(amt.txt-1),txt, adj=0)\n",
 203     "           }}\n",
 204     "       }\n",
 205     "       else print('no spikes or more than one spike')\n",
 206     "     }\n",
 207     "     invisible(list(heights = heights, xbr = xbr))\n",
 208     "   }\n",
 209     "   else {\n",
 210     "     return(list(heights = heights, xbr = xbr,counts=counts))\n",
 211     "   }\n",
 212     " }\n",
 213     "\n",
 214     "#' Calculate interquartile range\n",
 215     "#'\n",
 216     "#' Calculates the 25th and 75th quantiles given a vector x; used in function \\link{dhist}.\n",
 217     "#' @title Interquartile range\n",
 218     "#' @param x vector\n",
 219     "#' @return numeric vector of length 2, with the 25th and 75th quantiles of input vector x. \n",
 220     "iqr<-function(x){\n",
 221     "   return(diff(quantile(x, c(0.25, 0.75), na.rm = TRUE)))\n",
 222     "}"
 223    ]
 224   },
 225   {
 226    "cell_type": "markdown",
 227    "metadata": {},
 228    "source": [
 229     "## Reading dataframe"
 230    ]
 231   },
 232   {
 233    "cell_type": "code",
 234    "execution_count": 6,
 235    "metadata": {},
 236    "outputs": [
 237     {
 238      "name": "stdout",
 239      "output_type": "stream",
 240      "text": [
 241       "       X               op               msg_size     start        \n",
 242       " Min.   :  1130   Length:3500        Min.   :1   Min.   : 0.2003  \n",
 243       " 1st Qu.: 82565   Class :character   1st Qu.:2   1st Qu.:14.3932  \n",
 244       " Median :179094   Mode  :character   Median :4   Median :30.8789  \n",
 245       " Mean   :174237                      Mean   :4   Mean   :30.3998  \n",
 246       " 3rd Qu.:262094                      3rd Qu.:6   3rd Qu.:45.7296  \n",
 247       " Max.   :358579                      Max.   :7   Max.   :63.0828  \n",
 248       "    duration          experiment            type               index       \n",
 249       " Min.   :1.510e-07   Length:3500        Length:3500        Min.   :  2260  \n",
 250       " 1st Qu.:1.560e-07   Class :character   Class :character   1st Qu.:165130  \n",
 251       " Median :1.580e-07   Mode  :character   Mode  :character   Median :358189  \n",
 252       " Mean   :2.322e-07                                         Mean   :348475  \n",
 253       " 3rd Qu.:2.270e-07                                         3rd Qu.:524188  \n",
 254       " Max.   :1.452e-05                                         Max.   :717158  \n",
 255       "    origin         \n",
 256       " Length:3500       \n",
 257       " Class :character  \n",
 258       " Mode  :character  \n",
 259       "                   \n",
 260       "                   \n",
 261       "                   \n"
 262      ]
 263     },
 264     {
 265      "data": {
 266       "text/html": [
 267        "<table class=\"dataframe\">\n",
 268        "<caption>A data.frame: 6 × 9</caption>\n",
 269        "<thead>\n",
 270        "\t<tr><th></th><th scope=col>X</th><th scope=col>op</th><th scope=col>msg_size</th><th scope=col>start</th><th scope=col>duration</th><th scope=col>experiment</th><th scope=col>type</th><th scope=col>index</th><th scope=col>origin</th></tr>\n",
 271        "\t<tr><th></th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
 272        "</thead>\n",
 273        "<tbody>\n",
 274        "\t<tr><th scope=row>1</th><td>1130</td><td>MPI_Send</td><td>3</td><td>0.200334</td><td>2.57e-07</td><td>grenoble</td><td>exp/exp_PingPong.csv</td><td>2260</td><td>real</td></tr>\n",
 275        "\t<tr><th scope=row>2</th><td>1131</td><td>MPI_Send</td><td>3</td><td>0.200338</td><td>1.85e-07</td><td>grenoble</td><td>exp/exp_PingPong.csv</td><td>2262</td><td>real</td></tr>\n",
 276        "\t<tr><th scope=row>3</th><td>1132</td><td>MPI_Send</td><td>3</td><td>0.200342</td><td>1.63e-07</td><td>grenoble</td><td>exp/exp_PingPong.csv</td><td>2264</td><td>real</td></tr>\n",
 277        "\t<tr><th scope=row>4</th><td>1133</td><td>MPI_Send</td><td>3</td><td>0.200345</td><td>1.56e-07</td><td>grenoble</td><td>exp/exp_PingPong.csv</td><td>2266</td><td>real</td></tr>\n",
 278        "\t<tr><th scope=row>5</th><td>1134</td><td>MPI_Send</td><td>3</td><td>0.200349</td><td>1.54e-07</td><td>grenoble</td><td>exp/exp_PingPong.csv</td><td>2268</td><td>real</td></tr>\n",
 279        "\t<tr><th scope=row>6</th><td>1135</td><td>MPI_Send</td><td>3</td><td>0.200352</td><td>1.64e-07</td><td>grenoble</td><td>exp/exp_PingPong.csv</td><td>2270</td><td>real</td></tr>\n",
 280        "</tbody>\n",
 281        "</table>\n"
 282       ],
 283       "text/latex": [
 284        "A data.frame: 6 × 9\n",
 285        "\\begin{tabular}{r|lllllllll}\n",
 286        "  & X & op & msg\\_size & start & duration & experiment & type & index & origin\\\\\n",
 287        "  & <int> & <chr> & <int> & <dbl> & <dbl> & <chr> & <chr> & <int> & <chr>\\\\\n",
 288        "\\hline\n",
 289        "\t1 & 1130 & MPI\\_Send & 3 & 0.200334 & 2.57e-07 & grenoble & exp/exp\\_PingPong.csv & 2260 & real\\\\\n",
 290        "\t2 & 1131 & MPI\\_Send & 3 & 0.200338 & 1.85e-07 & grenoble & exp/exp\\_PingPong.csv & 2262 & real\\\\\n",
 291        "\t3 & 1132 & MPI\\_Send & 3 & 0.200342 & 1.63e-07 & grenoble & exp/exp\\_PingPong.csv & 2264 & real\\\\\n",
 292        "\t4 & 1133 & MPI\\_Send & 3 & 0.200345 & 1.56e-07 & grenoble & exp/exp\\_PingPong.csv & 2266 & real\\\\\n",
 293        "\t5 & 1134 & MPI\\_Send & 3 & 0.200349 & 1.54e-07 & grenoble & exp/exp\\_PingPong.csv & 2268 & real\\\\\n",
 294        "\t6 & 1135 & MPI\\_Send & 3 & 0.200352 & 1.64e-07 & grenoble & exp/exp\\_PingPong.csv & 2270 & real\\\\\n",
 295        "\\end{tabular}\n"
 296       ],
 297       "text/markdown": [
 298        "\n",
 299        "A data.frame: 6 × 9\n",
 300        "\n",
 301        "| <!--/--> | X &lt;int&gt; | op &lt;chr&gt; | msg_size &lt;int&gt; | start &lt;dbl&gt; | duration &lt;dbl&gt; | experiment &lt;chr&gt; | type &lt;chr&gt; | index &lt;int&gt; | origin &lt;chr&gt; |\n",
 302        "|---|---|---|---|---|---|---|---|---|---|\n",
 303        "| 1 | 1130 | MPI_Send | 3 | 0.200334 | 2.57e-07 | grenoble | exp/exp_PingPong.csv | 2260 | real |\n",
 304        "| 2 | 1131 | MPI_Send | 3 | 0.200338 | 1.85e-07 | grenoble | exp/exp_PingPong.csv | 2262 | real |\n",
 305        "| 3 | 1132 | MPI_Send | 3 | 0.200342 | 1.63e-07 | grenoble | exp/exp_PingPong.csv | 2264 | real |\n",
 306        "| 4 | 1133 | MPI_Send | 3 | 0.200345 | 1.56e-07 | grenoble | exp/exp_PingPong.csv | 2266 | real |\n",
 307        "| 5 | 1134 | MPI_Send | 3 | 0.200349 | 1.54e-07 | grenoble | exp/exp_PingPong.csv | 2268 | real |\n",
 308        "| 6 | 1135 | MPI_Send | 3 | 0.200352 | 1.64e-07 | grenoble | exp/exp_PingPong.csv | 2270 | real |\n",
 309        "\n"
 310       ],
 311       "text/plain": [
 312        "  X    op       msg_size start    duration experiment type                \n",
 313        "1 1130 MPI_Send 3        0.200334 2.57e-07 grenoble   exp/exp_PingPong.csv\n",
 314        "2 1131 MPI_Send 3        0.200338 1.85e-07 grenoble   exp/exp_PingPong.csv\n",
 315        "3 1132 MPI_Send 3        0.200342 1.63e-07 grenoble   exp/exp_PingPong.csv\n",
 316        "4 1133 MPI_Send 3        0.200345 1.56e-07 grenoble   exp/exp_PingPong.csv\n",
 317        "5 1134 MPI_Send 3        0.200349 1.54e-07 grenoble   exp/exp_PingPong.csv\n",
 318        "6 1135 MPI_Send 3        0.200352 1.64e-07 grenoble   exp/exp_PingPong.csv\n",
 319        "  index origin\n",
 320        "1 2260  real  \n",
 321        "2 2262  real  \n",
 322        "3 2264  real  \n",
 323        "4 2266  real  \n",
 324        "5 2268  real  \n",
 325        "6 2270  real  "
 326       ]
 327      },
 328      "metadata": {},
 329      "output_type": "display_data"
 330     }
 331    ],
 332    "source": [
 333     "df = read.csv(text=df_str)\n",
 334     "print(summary(df))\n",
 335     "head(df)"
 336    ]
 337   },
 338   {
 339    "cell_type": "markdown",
 340    "metadata": {},
 341    "source": [
 342     "## Clustering"
 343    ]
 344   },
 345   {
 346    "cell_type": "markdown",
 347    "metadata": {},
 348    "source": [
 349     "Calculating the intercept for the duration of each message."
 350    ]
 351   },
 352   {
 353    "cell_type": "code",
 354    "execution_count": 7,
 355    "metadata": {},
 356    "outputs": [],
 357    "source": [
 358     "df$intercept_estimate = df[[metric]] - df[[coeff_col_name]]*coeff_value"
 359    ]
 360   },
 361   {
 362    "cell_type": "markdown",
 363    "metadata": {},
 364    "source": [
 365     "Running dhist."
 366    ]
 367   },
 368   {
 369    "cell_type": "code",
 370    "execution_count": 55,
 371    "metadata": {},
 372    "outputs": [
 373     {
 374      "name": "stderr",
 375      "output_type": "stream",
 376      "text": [
 377       "Warning message in hist.default(x, breaks = cut.pt, plot = FALSE, probability = TRUE):\n",
 378       "“argument ‘probability’ is not made use of”\n"
 379      ]
 380     },
 381     {
 382      "data": {
 383       "text/plain": [
 384        "$heights\n",
 385        "\n",
 386        "          1           2           3           4           5           6 \n",
 387        "4510.135322 5863.669177 7666.843201 2453.222423 1045.499552  347.102832 \n",
 388        "          7           8           9          10          11          12 \n",
 389        "1448.764591  211.473846    0.000000    0.000000    0.000000    0.000000 \n",
 390        "         13 \n",
 391        "   1.778268 \n",
 392        "\n",
 393        "$xbr\n",
 394        " [1] 1.054988e-07 1.210706e-07 1.354500e-07 1.482141e-07 1.828565e-07\n",
 395        " [6] 2.538565e-07 4.000706e-07 5.307136e-07 8.678415e-07 1.523917e-06\n",
 396        "[11] 2.675975e-06 4.698973e-06 8.251327e-06 1.447928e-05\n",
 397        "\n",
 398        "$counts\n",
 399        "\n",
 400        "       1        2        3        4        5        6        7        8 \n",
 401        "620.9278 658.0722 690.4405 515.2837 342.9944 157.8873 409.3941 104.0000 \n",
 402        "       9       10       11       12       13 \n",
 403        "  0.0000   0.0000   0.0000   0.0000   1.0000 \n"
 404       ]
 405      },
 406      "metadata": {},
 407      "output_type": "display_data"
 408     }
 409    ],
 410    "source": [
 411     "data = df$intercept_estimate\n",
 412     "if (log) {\n",
 413     "    data = log(data)\n",
 414     "}\n",
 415     "\n",
 416     "dh = dhist(data, plot=FALSE)\n",
 417     "\n",
 418     "dh"
 419    ]
 420   },
 421   {
 422    "cell_type": "markdown",
 423    "metadata": {},
 424    "source": [
 425     "Saving result to cvs"
 426    ]
 427   },
 428   {
 429    "cell_type": "code",
 430    "execution_count": 56,
 431    "metadata": {},
 432    "outputs": [],
 433    "source": [
 434     "write.csv(list(xbr=dh$xbr), xbr_file)\n",
 435     "write.csv(dh$heights, height_file)"
 436    ]
 437   },
 438   {
 439    "cell_type": "markdown",
 440    "metadata": {},
 441    "source": [
 442     "## Visualizing"
 443    ]
 444   },
 445   {
 446    "cell_type": "markdown",
 447    "metadata": {},
 448    "source": [
 449     "Histogram with vertical lines at the mean of each cluster."
 450    ]
 451   },
 452   {
 453    "cell_type": "code",
 454    "execution_count": 58,
 455    "metadata": {},
 456    "outputs": [
 457     {
 458      "data": {
 459       "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAC4lBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8RERETExMUFBQVFRUWFhYXFxcY\nGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIkJCQmJiYnJycoKCgpKSkqKiosLCwt\nLS0uLi4vLy8xMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9A\nQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NU\nVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVm\nZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4\neHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OFhYWGhoaHh4eIiIiJiYmKioqL\ni4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2e\nnp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+w\nsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHC\nwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU\n1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm\n5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4\n+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///9F8m1nAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElE\nQVR4nO3de5xcZYHm8aMMKmBYHTXLis6oeGNixFEjF8UBxySgIgIaSBBXgdEM4OKKIgyEGVBR\nBy8RFNYri3gZgmZQENFFVhxlFHWQy6AEMVS6w8WI4ZL6f6u6uvvU83S68lRz+nRV9vf9fKhT\np85b5+23c37pk1tTNAE8asVcfwDA9oCQgAoQElABQgIqQEhABQgJqAAhARUgJKACFYX03iW9\nLH7t4iVL9n3hC/Zc+KIlr1zw6vZL+y34m9bjgQv2kYGvWbB3zxO9aK+eh91rxzcLXtbX2/r2\n0gWLu/bGVjsHXrvtIfssOLDqWQd4tY/CtBearfaYqkN62297Hd3UeLDZ/EBRFH+xU/MjxRfb\nL51VfK31eHVxsgz8SfGOntPs+vy+PqrRzqZRvLavt/Xt1cX9XXsPNB6Y3emmMbLtIScWP6h6\n1gcbm6o+ZeSeh2fz7NNdaA837uvefXDJ5FNCqgAh1Y+QCGmWEFJ1CEkQkiGkECEJQjKEFCIk\nQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRB\nSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFI\nhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiG\nkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUGOaQXtELaZYfDCGm2EFJ1CEkQ\nkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCS\nIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIh\npBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGk\nECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQ\nIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAh\nCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGF5iykY/7z4R7+\n0PjTww+Ph/Th4vPtl84sLm09XlX8Dxn44+Ltvc7z8K7P73nYjXY2vy/+tq+39W3/4p6uvU2N\nTbM73TRGtj3khOKaqmf9U+MPVZ8ysnHzbJ59ugttc+Pe7t0HKg/p6Bs39jDSGN24cTyks4rz\n2y99oPh863FNsVIGfq9Y0es8G+ft0fOw29DZ3FIc2Nfb+rZfsa5rb6QxMrvTTaOx7SHHF2ur\nnnV0jla7YXQ2zz7dhTba2NC9e3d5q8OtXQW4tavfdnprR0glQqoDIRHSLCGk6hCSICRDSCFC\nEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIS\nhGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKE\nZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRk\nCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQI\nKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgp\nREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClE\nSIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURI\ngpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiC\nkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQ\nDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAM\nIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwh\nhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGF\nCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUI\nSRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJ\nEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFqg9p\n5Jw3v/E9NzWb93/oyMNPXz91O46QSoRUh2EL6e9PvvV3HzzigeYZJ9925znHPTJlO46QSoRU\nhyEL6b5VrULuXvzrxpJbW1+FDrrBtxPjCKlESHUYspDG/Grp6LWv39J6cvzFvp0YQkglQqrD\nEIZ037Gfba49qv3slNW+bT1ctbDljf/e2JbxkM4sPtXee19xUevxG8U7ZcxVxfKe55i3xzan\n2YqbigNm8rbcvsXtsztBVY4tLp/rD2E4hBfa78qfobcd0h1v+/iW5trl4wHZtvVwwztaDv/F\nPT2MNjbec894SP9YnN9+6QPFF1qPlxcrZeA1xdG9znPPvOf2POwanc2txWv6elvf9ivu7Nob\nbYzO7nTTaGx7yPHF2qpn3ThHq92wcTbPPt2FtrEx0r3b6COkGw6/rPV4XedW7hLfTozi1q7E\nrV0dhu3W7heH/bi9GVlyc7N579IbfTsxjJBKhFSHIQtp8zFfat8LPtA8a+Vt6047YcuU7ThC\nKhFSHYYspBsWj1nT3HTusiNWtS5M344jpBIh1WHIQkoRUomQ6kBIhDRLCKk6hCQIyRBSiJAE\nIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQh\nGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZ\nQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlC\nChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIK\nEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoR\nkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGS\nICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIg\nJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAk\nQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRD\nSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENI\nIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0gh\nQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFC\nEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIS\nhGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKE\nZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSKG5C+k3\nW3rY1Ni8Zct4SOcWX2i/tKr4auvxu8XJMvDfinf0Os+WXZ/f87Ab7WzuLl7b19v69urivq69\nPzb+OLvTTWNk20NOLL5f9aybG5uqPmXknodm8+zTXWgPNbp/pLdsrjyko38+2sNIY2R0dDyk\nVcWn2y+dWnyu9XhZ8S4ZeHWxotd5Ruft0fOw29DZ3Fwc2Nfb+rZvcUfX3thq50Bj20OOK75V\n/bRzs9oNszrtdBfaSGND9+768laHW7sKcGtXv+311o6QJhFSHQiJkGYJIVWHkAQhGUIKEZIg\nJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAk\nQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBCjzKkhb/sbL/yvH4nJqQSIdVh\noEMqrh/bPHT64/qdmJBKhFSHAQ6pKL2o34kJqURIdRjgkG74aLH06La3nnpHvxMTUomQ6jDA\nITWbB/56phMTUomQ6jDQIc0cIZUIqQ4DHdL6I3d7bOcXSf1OTEglQqrDQId0yJ/tf+TYr5KO\n7ndiQioRUh0GOqQnf32mExNSiZDqMNAh7Xz3TCcmpBIh1WGgQ9r7uzOdmJBKhFSHgQ7pxy+5\ndoYTE1KJkOow0CEt2r3Y+Zlj+p2YkEqEVIeBDmnv/Sf0OzEhlQipDgMd0swRUomQ6kBIhDRL\nCKk6j/bPkSbM63diQioRUh0GOqSlY16y0wuO63diQioRUh0GOqRxd+2zpt+JCalESHUYhpCa\n1y/sd2JCKhFSHYYipLt26ndiQioRUh2GIaQtZz6934kJqURIdRjokP5qzAv+vDip34kJqURI\ndRiCkBa86qOb+52YkEqEVIeBDmnmCKlESHUY8JA2rFl9/tr7mn0jpBIh1WGgQ3rkxB3b37Bh\nl7P7npiQSoRUh4EO6ezi4Au+ueZTBxYX9TsxIZUIqQ4DHdLzTuhs38Z3Wu0LIdVvoEN6/JWd\n7eX8gWxfCKl+Ax3SLpd1tl9/Yr8TE1KJkOow0CG94pVjf4D0wAH79TsxIZUIqQ4DHdLlj3nG\n28/4h2N2e+y3+52YkEqEVIeBDqn5tee2f/v7hZf3PTEhlQipDoMdUrN554+u//0MJiakEiHV\nYbBDuutjrYe7T1/f98SEVCKkOgx0SP8xv/3/vLy9mH9rvxMTUomQ6jDQIR307B+1N7989uv6\nnZiQSoRUh4EO6Smf6Ww/xXcR6gsh1W+gQ9rp853tF3bud2JCKhFSHQY6pJcfOPbR3bfXon4n\nJqQSIdVhoENa+5hnHXfaqcuf8ti1/U5MSCVCqsNAh9S8YmH7D2T35A9k+0NI9RvskJrNDT/7\nxQz+gSwhde0RUh0GPaQZIqQSIdWBkAhplhBSdQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBS\niJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKI\nkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQ\nBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAE\nIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQh\nGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZ\nQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlC\nChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIK\nEZIgJENIoVkIad2JS9ub+z905OGnr5+6HUdIJUKqw7CFdM2yc8dCOuPk2+4857hHpmzHEVKJ\nkOowbCFdefcP2yE1ltza+ip00A2+nRhGSCVCqsOwhdRsjoV07eu3tB6Pv9i3E4MIqURIdRjS\nkNYe1X56ymrfth6+/8qWQ3+2oZdG67/xkM4sVrdfeX/xv1qP/1K8U8Z9t1je8zzz9uh5eGvT\nttxUHNDX2/q2b/EbmbYxu9NNI5j12OLyuZh2NszutNNeaPpje1f5M3Qa0vLxgGzberj+zS2H\n3bixh5HG6MaN4yGdVZzffukDxedbj2uKlTLwe8WKXufZOG+Pnofdhs7mluLAvt7Wt/2KdV17\nI42R2Z1uGo1tDzm+WFv1rKNztNoNo7N59ukutNHGhu7du/sO6brOrdwlvp0YxK1diVu7Ogzp\nrd3IkpubzXuX3ujbiUGEVCKkOgxbSKONK5Y2WlfGWStvW3faCVumbMcRUomQ6jBsIR29uO0b\nzU3nLjtiVevC9O04QioRUh2GLaQQIZUIqQ6EREizhJCqQ0iCkAwhhQhJEJIhpBAhCUIyhBQi\nJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIk\nQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRB\nSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFI\nhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiG\nkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQ\nQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBC\nhCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKE\nJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQk\nCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQI\nyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJ\nEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQ\nUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBS\niJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKI\nkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQ\nBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBCcxbSip9t6KXR\n+m88pDOL1e1X/qrY+7DDDtu/eKeM+26xvOd55u3R8/A0bioOmMnbcvsWv+nebTRmd7ppBLMe\nW1w+F9POhtn9JE97oem0d5U/Q8/VV6SJkPiKVBm+IlVneG7tCKlyhFQdQhKEZAgpREiCkAwh\nhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGF\nCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUI\nSRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJ\nEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQ\nkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCS\nIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIh\npBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGk\nECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQ\nIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAh\nCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSKEhCOlFxctbBR1GSJUjpOoQkiAkQ0gh\nQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFC\nEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIS\nhGQIKURIgpAMIYUISRCSIaQQIQlCMoQUIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKE\nZAgpREiCkAwhhYYvpOe3dyYHEtKjQEjVISRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKE\nZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQiJEFIhpBChCQIyRBSaEhDGkdIjwohVYeQBCEZQgoR\nkiAkQ0ghQhKEZAgpREiCkAwhheoL6f4PHXn46esndwmpREh12F5COuPk2+4857hHJnYJqURI\nddhOQmosubX1VemgGyb2CalESHXYTkK69vVbWo/HXzyxT0glQqrDdhLS2qPaj6esbj38cEnL\nYT8f7WGkMTI6uu2QRkevLlb0Os/ovD16HnYbOpubiwP7elvf9i3u6NobW+0c2LDtIccV36p8\n2sFd7aMw3YU20pB515c/Q880pOV9h3RyK6RnPGF0VfHp9kunFp9rPV5WvEsGEtKjQEjVqS2k\n6zq3dpdM7Ae3dh9ohfQXOzU/Unyx/dJZxddaj1fz7biqw61ddWq7tRtZcnOzee/SGyf2CalE\nSHXYTkJqnrXytnWnnbBlYpeQSoRUh+0lpE3nLjti1ejkLiGVCKkO20tIhpBKhFQHQiKkWUJI\n1SEkQUiGkEKEJAjJEFKIkAQhGUIKEZIgJENIIUIShGQIKURIgpAMIYUISRCSIaQQIQlCMoQU\nIiRBSIaQQoQkCMkQUoiQBCEZQgoRkiAkQ0ghQhKEZAgpREiCkAwhhQhJEJIhpBAhCUIyhBQi\nJEFIhpBChCQIyRBSiJAEIRlCChGSICRDSCFCEoRkCClESIKQDCGFCEkQkiGkECEJQjKEFCIk\nQUiGkEKEJAjJEFKIkAQhGUIKzV1IH72wh0+fd8GFF755/vz5z9v9wmPm/337paPmv7v1eOr8\nQ2XgqvkH9TrPhc/cs+dh94nO5rz5r+zrbX17xfzVXXvnn3f+7E43jY9ve8gh80+vetYLzvt0\n1aeMfPIzs3n26S60z5z3qe7dCyoP6XuX9nLmitWXXvrh173ikHeecOk/r/xk+6WPrGy9dOln\nVn5QBl608p96nujE9/Y87C7ubL688h/6elvfTlv5v7v2Przi3NmdbhoXb3vIOSs/W/Wsn1xx\nVtWnjFz8ldk8+3QX2kUr3if7/1p1SL2dt/BHdUwzIC5d+I25/hBqdN3CT8z1h1Cj3yx8/zRH\nCKlyhLT9IqQaEdL2i5BqREjbrzkOCdjeERJQAUICKkBIQAUqDen+Dx15+Onrpz7f6qBpBwyP\nv1vcckjn+VaXM/ni5W89+Pih/w2XdScunXjac7Xdn5bh1bXarS7IVltpSGecfNud5xz3yJTn\n4zPfJAemDhg6yy9rNBrjf8Vtq8uZePE7y65f//Vj5uYvpVXmmmXnTl5aPVfb/WkZWt2r3eqC\nbLVVhtRYcmsrl4NukOejZy97w3tuab/0s3d2H+gePKzecP3k062st+vFY66co4+wSlfe/cOJ\nS6v3ars+LcOra7WTC+q12ipDuvb1W1qPx18sz088+77Nn3vz5uZESBMHugcPqQcXf+xdK1at\nG3u+lfWWL25YfOXfveHEX83dR1qRyUur52q7Py3DbHK1kwvqtdoqQ1p7VPvxlNXdz29ZPNps\nbnnTNc2JkCYOdA8eUve85cM33XTaW/7Qfr6V9ZYv3rT4f95x3+o33TNnH2lFJi+tnqvt/rQM\ns8nVTiyo52orDWl559zdz69ZPOaSGw499JAlhx56wuSB7sHD7I+HXNHe+Hq/v3Tp0l9OvHjT\n4tYd7MOHfWcuP84qlCH1Wu3YiPFPyzArb+3aWgvqudoqQ7qu89Xuku7n1y3ufCXcvH79Nceu\nX7+hPNA1eKgdO/YPFX29m26//fY/TbzYWHxza3vc0K918tLqudrOkM6nZZhpSK0F9VxtlSGN\nLGldL/cuvbH7+W8X/0dr7672S51bu4kD3YOH1O3//FCz+cAhV7Wfb2W95YuPLLus9VPJG6+Z\nw4+1EpOXVs/Vdn9ahtnkaicW1HO1lf7291krb1t32glbmlf8S/n8lGvZR5IAAAR2SURBVHff\n/fA339D+vcNOSJMHJrbD677Dz71r3arlf5pmvc3yxUuO+GnjY8vm5l+gV2a0ccXSRuOBba52\n8tMy1LpWO7mgXqutNKRN5y47YlXrF2Rnv698PvpPh77x3TduZdDEdojd+r5D33zG73usd+LF\nRy56y8Hv6fltLYbA0WO/QvjGtlc78WkZat2rnVhQr9XyV4SAChASUAFCAipASEAFCAmoACEB\nFSAkoAKEBFSAkObMX++xPU3z/ztCmjPnrure+2nfPxLbfsfYCJ2m73MgwydyQHys7x+Jbb+j\nihHI8ImcM+17rr1f8ZNXzXvKm9Y3DyyKYmGzefWr5+204ILWwUV7X/b0lzWbV+zzxKcd0v43\nGJMHXvTSK/fa6UnL75l4R2lyyO/e+ozHP+11vxof0Znmmr2esNvZD5682xP3v7U14kt77TRv\n4ZeaU2bFjBHSnGlf4fvvvte3139lhyObv15aXP/L5nd22OeyK95efLDZfNWezz1vTfOKxxzw\n+Que9V/v6jrwsqe8+P80PrfjwePvKJVDXjr//Ku+8MKnbuqMGJvm6a/8tzsOLl59+rrv7fq3\nzeaXi4PXrHlNsWbKrJgxQpozY1f42P9Cb//dms2j2z8SC57d/k5DS+Y90Drw1dazF//lQ83m\n/33cR7sOLCra/6zp6OK3nXeUJofcW7yn9eSWVXd2RnSmuaHZ/H7x8tbuEbs0m6tetbnZvPfP\njpgyK2aMkObM2BW+c/vZkY/tXNLri3c90PLJ4kfN/R/3YLO5oTi2M7TrwKJd2v+E66LimxZS\nOeTBJz/zO51vlFWG1KqneUtxUuvxpGLif9749L2nzIoZI6Q5M3aFP7P9rH05t//7aTHuq2Nf\npJo/L07rDO06sOhZ7RfWFJ+1kLqG/OAviye//gsPdYfUnuY/i7NajycXG5v3vv8Fu+6wQ7Fo\nyqyYMUKaM1sLacUPxzQ6B35RnNoZ2nWgE9LXi4umhDQ5pPnwlSc9v3jxH6cNaZ8d3nvNz36+\n26Ips2LGCGnOTA1ppDhy4uDYgfuKsW9Vc/vdXQcW7dT+33l/qlhrIXUNGfPx4sLpQrq5OKb1\n5KEnLJoyK2aMkOaMhvTWonUv9pL/srG1e9EpD40feOFTWr+g+VXrBq88sKj1q6Nm86DHj3be\nUZoc8uND29+U+5binM6IqSH9sji92f4jpJdOmRUzRkhzRkM6tTj9K82rd9zzon99345HTRxY\n89i//uLq5zz1rq4Di3Z/zie+fVKxbPwdpckhd83b84Jvf/nlu97SGTE1pAd3/2/f+MGJ++03\n76o/2KyYMUKaMxrSHQt2bO1//2/m7ficsx+aONC8/KU7P/XgXze7Dix67o/32flJx9w/8Y7S\n5JB/P/ipO+528E/GR2zl10jXv2znp/33ey/78yfdZLNixghpyCzi76AOJEIaMoQ0mAhpyBDS\nYCKkISMhfauY9Ik5+4jQRkjD7P6fTxru71k7/AgJqAAhARUgJKAChARUgJCAChASUAFCAirw\n/wBOznRqSUaAyQAAAABJRU5ErkJggg==",
 460       "text/plain": [
 461        "plot without title"
 462       ]
 463      },
 464      "metadata": {
 465       "image/png": {
 466        "height": 420,
 467        "width": 420
 468       }
 469      },
 470      "output_type": "display_data"
 471     }
 472    ],
 473    "source": [
 474     "library(ggplot2)\n",
 475     "\n",
 476     "if (log) {\n",
 477     "    dh$xbr = exp(dh$xbr)\n",
 478     "}\n",
 479     "\n",
 480     "ggplot(df, aes(x=intercept_estimate)) + geom_histogram(bins=100) +\n",
 481     "    geom_vline(xintercept=dh$xbr) + \n",
 482     "    theme_bw()"
 483    ]
 484   },
 485   {
 486    "cell_type": "markdown",
 487    "metadata": {},
 488    "source": [
 489     "Same but in log scale"
 490    ]
 491   },
 492   {
 493    "cell_type": "code",
 494    "execution_count": 59,
 495    "metadata": {},
 496    "outputs": [
 497     {
 498      "name": "stderr",
 499      "output_type": "stream",
 500      "text": [
 501       "Warning message:\n",
 502       "“Transformation introduced infinite values in continuous y-axis”\n",
 503       "Warning message:\n",
 504       "“Removed 59 rows containing missing values (geom_bar).”\n"
 505      ]
 506     },
 507     {
 508      "data": {
 509       "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAC6FBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8RERESEhITExMUFBQVFRUWFhYX\nFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIkJCQmJiYnJycoKCgpKSkqKior\nKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9\nPT0+Pj5AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09Q\nUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFi\nYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0\ndHR1dXV3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn6AgICBgYGCgoKDg4OFhYWGhoaHh4eIiIiJ\niYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqb\nm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKyt\nra2urq6vr6+wsLCxsbGzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/A\nwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS\n0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk\n5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX2\n9vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///+V7mQAAAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nO3df5xdZX3g8SNptEkErUWzKNoqaxakjSFITMOvWkKbSUChGJHwI2Kh2DUL\nFCssIAsWatOgLCTCKlVbLWxXCDYNxrbBgjW2soBI+LUshI3JzQwCYkogOf/ufb73uc+cZ77n\n3nme771xMjOfz+vl3HnOOfecb27O27nJJKEoiajnirEegGgiBCSiPgQkoj4EJKI+BCSiPgQk\noj4EJKI+BCSiPtQjpE8v7tKiRfLwu7Pnq13vmx0tB2YfNbw4fvZxXU46973dLjmyebMHuu0e\n6Lq32ZGzRzmg1aLZ81In6tJ7j6w58yLz6Y7MeqVGbd7s36vdbhxwweyFPQyjOnb28Z12DSwa\n7WdZd3zn00Wd2y9IH3+6y84XGrvdw6PFGWrX3CnRclexYHjxl8UXupx01hty5vtA8UK33Ttf\nGuX5c6YmXWZ7MZA6UZemH6427Wn81Hy6w17fyzCqxcW22u3P7TadbkVxXy/TjOz64quddu1s\n/Hv26b5YfDHlsF2Lw6dA6hqQQkDSAckHpPSApAOSD0jpAUkHJB+Q0gOSDkg+IKUHJB2QfEBK\nD0g6IPmAlB6QdEDyASk9IOmA5ANSekDSAckHpPSApAOSD0jpAUkHJB+Q0gOSDkg+IKUHJB2Q\nfEBKD0g6IPmAlB6QdEDyASk9IOmA5ANSekDSAckHpPSApAOSD0jpAUkHJB+Q0gOSDkg+IKUH\nJB2QfEBKD0g6IPmAlB6QdEDyASk9IOmA5ANSekDSAckHpPSApAOSD0jpAUkHJB+Q0gOSDkg+\nIKUHJB2QfEBKD0g6IPmAlB6QdEDyASk9IOmA5ANSekDSAckHpPSApNtHIS1tBaThgJQYkFoD\nAakEUi8BqTUQkEog9RKQWgMBqQRSLwGpNRCQSiD1EpBaAwGpBFIvAak1EJBKIPUSkFoDAakE\nUi8BqTUQkEog9RKQWgMBqQRSLwGpNRCQSiD1EpBaAwGpBFIvAak1EJBKIPUSkFoDAakEUi8B\nqTUQkEog9RKQWgMBqQRSLwGpNRCQSiD1EpBaAwGpBFIvAak1EJBKIPUSkFoDAakEUi8BqTUQ\nkEog9RKQWgMBqQRSLwGpNRCQSiD1EpBaAwGpBFIvAak1EJBKIPUSkFoDAakEUi8BqTUQkEog\n9RKQWgMBqQRSLwGpNVAV0rziCPfQ4aRASg1IrcYHpP+7p3MvNF51D5uLM9SuuVOi5cvFgubH\nGFKHk856Q5crqj5QPN9t989/Nsrz50xNusy2YiB1oi5NP1xt2t34qfl0h72+l2FUi4uf1G5/\n7lXT6VYU9/YyzchWFV/ptOvnjZ3Zp1tTrEk57OW+QVr+4FDndjTkYVNxmto1e0q03FbMa36M\nIXU46SEHdLmi6pji6W67B3eM8vzfnJp0mUeLhakTdWnaoXpbY7QJOzdrRg+z6E4sNtdub9hO\nd36xvodhVJ8tVnfaNdgYzD7dqmJVymHbFvULUn/e2i09rThwaYi3dsPx1i6xcf/WDkguIGUH\npCggSUDKDkhRQJKAlB2QooAkASk7IEUBSQJSdkCKApIEpOyAFAUkCUjZASkKSBKQsgNSFJAk\nIGUHpCggSUDKDkhRQJKAlB2QooAkASk7IEUBSQJSdkCK6hFS2w6QXECyN1kh/cprllYDkgtI\n9oAEpBCQ7AEJSCEg2QMSkEJAsjfpIHksQKoJSPaABKQQkOwBCUghINkDEpBCQLIHJCCFgGQP\nSEAKAckekIAUApI9IAEpBCR7QAJSCEj2gASkEJDsAQlIISDZAxKQQkCyByQghYBkD0hACgHJ\nHpCAFAKSPSABKQQke0ACUghI9oAEpBCQ7AEJSCEg2QMSkEJAsgckIIWAZA9IQAoByR6QgBQC\nkj0g7duQfmW/pV1mCQEpOyBFAUkCUnZAigKSBKTsgBQFJAlI2QEpCkgSkLIDUhSQJCBlB6So\ndEgRHCDVBSR7QAJSCEj2gASkEJDsAQlIISDZAxKQQkCyByQghYBkD0hACgHJHpCAFAKSPSAB\nKQQke0ACUghI9oAEpBCQ7AEJSCEg2QMSkEJAsgckIIWAZA9IQAoByR6QgBQCkj0gASkEJHtA\nAlIISPaABKQQkOwBCUghINkDEpBCQLIHJCCFgGQPSEAKAckekIAUApI9IAEpBCR7QAJSCEj2\ngASkEJDsAQlIISDZAxKQQkCyByQghYBkD0hACgHJHpCAFAKSPSABKQQke0ACUghI9oAEpBCQ\n7AEJSCEg2QPSuIDk63gYkLIDUhSQJCBlB6QoIElAyg5IUUCSgJQdkKKAJAEpOyBFAUkCUnZA\nigKSBKTsgBQFJAlI2QEpCkgSkLIDUhSQJCBlB6QoIElAyg5IUUCSgJQdkKKAJAEpOyBFAUkC\nUnZAigKSBKTsgBQFJAlI2QEpCkgSkLIDUhSQJCBlB6SoBEhLFxW/tnRkQKoJSPaABKQQkOwB\nCUghINkDEpBCQLIHJCCFgGQPSEAKAckekIAUApI9IAEpBCR7QAJSCEj2gASkEJDsAQlIISDZ\nAxKQQkCyByQghYBkD0hACgHJHpCAFAKSPSABKQQke0ACUghI9oAEpBCQ7AEJSCEg2QMSkEJA\nsgckIIWAZA9IQAoByR6QgBQCkj0gASkEJHtAAlIISPaABKQQkOwBCUghINkDEpBCQLK370Pa\nctES9/Di5878yJXbhh99QJKAlN1kg7Rx2UqBdNUlTz573QW7w6MPSBKQsptskDZsv89Baix+\novnV6KT724/t3UCSgJTdZINUlgLp3g/taX78xDfaj+2dQJKAlN0khbTuLPfppWvaj80P9xzb\n7LQHdnSp4T5YIXU45yEHdLviyI4pnuo6YGOU58eQOh62uTghZ6oOTTtUbxt1ws7NmtHDLLqF\nxSO1240Dnles62EY1TXFTR33GSZcWaxMOWzromxIZ3tIZwdImz7abOlDz3VusOE+WiF1OOkh\nB3S5ourY4pluuwcHR3l+DKnjYY8VC3Om6tC0Q/W2xg7z6WbN6GEW3YnFo7XbjQOeX9zdwzCq\nPy3WdNo11BjKPt2qYlXKYduzIX2v9ZbutvZjeydv7STe2mU3Sd/aDS5+rCyfX/JQ+7G9E0gS\nkLKbbJCGGuuXNBo7y2tWPLnligv3hEcfkCQgZTfZIC0fcN1RvrRy2elXD5Xh0QckCUjZTTZI\nowQkCUjZASkKSBKQsgNSFJAkIGUHpCggSUDKDkhRQJKAlB2QovYipA43L5BSA1IrIAEJSKkB\nCUghINkDEpBCQLIHJCCFgGQPSEAKAckekIAUApI9IAEpBCR7QAJSCEj2gASkEJDsAQlIISDZ\nAxKQQkCyByQghYBkD0hACgHJHpCAFAKSPSABKQQke0DKg7R06QFTu93RIwNS3wKSDkg+IKUH\nJB2QfEBKD0i68QQpvoeBlB6QWgEJSEBKDUhACgHJHpCAFAKSPSABKQQke0ACUghI9oAEpBCQ\n7AEJSCEg2QMSkEJAsgckIIWAZA9IQAoByR6QgBQCkj0gASkEJHtAAlIISPaABKQQkOwBCUgh\nINkDEpBCQLIHJCCFgGQPSEAKAckekIAUApI9IAEpBCR7QAJSCEj2gASkEJDsASkRUkggjXZj\nh8yQ2nMCqR2QdEDyASk9IOmA5ANSekDSAckHpPSApAOSD0jpAUkHJB+Q0gOSDkg+IKUHJB2Q\nfEBKD0i6iQJpdE9A6ltA0gHJB6T0gKQDkg9I6QFJByQfkNIDkg5IPiClByQdkHxASg9IOiD5\ngJQekHRA8gEpPSDpgOQDUnpA0gHJB6T0gKQDkg9I6QFJByQfkNIDkg5IPiClByQdkHxASg9I\nOiD5gJQekHRA8gEpPSDpgOQDUnpA0gHJB6T0gKQDkg9I6QFJByQfkNIDkg5IPiClByQdkHxA\nSg9IOiD5gJQekHRA8gEpPSDpJjSkaCOQ+haQdEDyASk9IOmA5ANSekDSAckHpPSApAOSD0jp\nAUk3MSHVbQRS3wKSDkg+IKUHJB2QfEBKD0g6IPmAlB6QdEDyASk9IOmA5ANSekDSAckHpPSA\npAOSD0jpAUk3wSDV1hoBSH0LSDog+YCUHpB0QPIBKT0g6YDkA1J6QNIByQek9ICkA5IPSOkB\nSQckH5DSA5IOSD4gpQckHZB8QEoPSDog+YCUHpB0QPIBKT0g6YDkA1J6QNIByQek9ICkA5IP\nSOkBSQckH5DSA5IOSD4gpQckHZB8QEoPSDog+YCUHpB0QPIBKT0g6YDkA1J6QNIByQek9ICk\nA5IPSOkBSQckH5DSA5JuMkBqNRNI/QpIOiD5gJQekHRA8gEpPSDpgOQDUnpA0gHJB6T0gKQD\nkg9I6QFJByQfkNIDkg5IPiClByQdkHxASg9IOiD5gJQekHRA8gEpPSDpgOQDUnpA0gHJB6T0\ngKQDkg9I6QFJByQfkNIDkq5/kJY/ONS5HQ33cYwhPd1lwKHBHZ321ELqeJ5Hi4XdLpPYtEP1\ntkbHCUdt1oweZtGdWGyu3d6wne78Yn0Pw6g+W6zutGuwMZh9ulXFqpTDti3qFyS+Ikl8RcqO\nr0hRQJKAlB2QooAkASk7IEUBSQJSdkCKApIEpOyAFAUkCUjZASkKSBKQsgNSFJAkIGUHpCgg\nSUDKDkhRQJKAlB2QooAkASk7IEUBSQJSdkCKApIEpOyAFAUkCUjZASkKSBKQsgNSFJAkIGUH\npCggSUDKDkhRQJKAlB2QooAkASk7IEUBSQJSdkCKApIEpOyAFAUkCUjZASkKSBKQsgNSFJAk\nIGUHpCggSUDKDkhRQJKAlB2QooAkASk7IEUBSQJSdkCKApIEpOyAFAUkCUjZTWhIcx5uPd4+\nK/WCQJKAlN2EhlRskodXrnxt6gWBJAEpuwkMqRjuvakXnMyQKhubkEazlhCQ7O1LkO5fVSxZ\n7vrY5c+kXhBIQAKSK3prd8KjuRcEEpCA5OJ37eoD0siApOsEaduZB+3X+kVS6gWBBCQguSIy\np/7S8WfKr5KWp14QSEACkiuC9KZv5l4QSEACkiuCNH177gWBBCQguSJIC/4h94JAAhKQXBGk\nHxx5b+YFgQQkILkiSPMPLqa/Q0q94L4P6ZRu93cmpBhLZQWk7CY0pAXHt0u9IJCABCTXhP+G\nLJD6FJB0QPIBKT0g6Tp+H6nd/qkXBBKQgOSKIC2Rjpx22AWpFwQSkIDkqntrt/XotakXBFIK\npG77KgHJ3j4Jqdw0J/WCQAISkFy1kLZOS70gkIAEJFcdpD3/7W2pFwQSkIDkiiD9hnTYrxYX\np14QSEACkqsG0uzjVr2cekEgAQlILr4hW197TiC1A5KuM6Qda9fcvK7rH5iOAxKQgOSKIO2+\naKr7BxtmXJt8QSABCUiuCNK1xcm3fGvt6hOKW1MvCCQgAckVQZp1Yevx4xPoX1oFUp8Ckq4T\npNdtaD3exTdk23MCqR2QdJ0gzbiz9fjN5J8CIAEJSK4I0m8dK99A2vk7x6ReEEhAApIrgnTX\na95+3lWfOfeg/e5OvSCQgAQkV/x9pP/1bvfb3++5K/mCQAISkFwj/2TDs9/f9JOMCwIJSEBy\nxZC2Xt/8sP3K+lesLiABCUiuCNIjM91/8/KpYuYTqReclJDi3NFAym5CQzrpXd93Dw+/64Op\nFwQSkIDkiiAd+D9aj6v5V4TacwKpHZB0nSBN88N8bXrqBYEEJCC5IkjvP+FV9/DC3PmpFwQS\nkIDkiiCte807L7ji8rMP3G9d6gWBBCQgueLf/l4/x31D9nC+IQskFZB0Xb4hu+OBH2X8BVkg\nAQlIEv9mQ33tOYHUDkg6IPmAlB6QdEDyASk9IOmA5ANSekDSAckHpPSApAOSD0jpAUkHJN/e\nhqSP7BaQ7AFJApILSPaAJAHJBSR7QJKA5AKSPSBJQHIByR6QJCC5gGQPSBKQXECyByRpPENy\nfbB4q/+s9iy1+3RAsgckCUguINkDkgQkF5DsAUkCkgtI9oAkAckFJHtAkoDkApI9IElAcgHJ\nHpAkILmAZA9I0sSB5IvPEm/sGJDsAUkCkgtI9oAkAckFJHtAkoDkApI9IElAcgHJHpAkILmA\nZA9IEpBcQLIHJAlILiDZA5IEJBeQ7AFJ2uuQOtzKQEoPSDog+YCUHpB0QPIBKT0g6YDkA1J6\nQNIByQek9ICkA5IPSOkBSQckH5DSA5IOSD4gpQckHZB8QEoPSDog+YCUHpB0QPIBKT0g6YDk\nA1J6QNIByQek9ICkA5IPSOkBSQckH5DSA5IOSD4gpQckHZB8QEoPSDog+YCUHpB0QPIBKT0g\n6YDkA1J6QNIByQek9ICkA5IPSOkBSQckH5DSA5IOSD4gpQckHZB8QEoPSDog+YCUHpB0Fkhb\nLlriHl783JkfuXL4FQWSBKTsJiukjctWCqSrLnny2esuCK8dkCQgZTdZIW3Yfp+D1Fj8RPOr\n0kn3tzcDSQJSdpMVUlkKpHs/tKf58RPfaG8EkgSk7CY5pHVnuU8vXdP8cM+xzU57YEeXGu7D\nvgRJDdjoNHp7TiOk+Cwdrj6yaYfWvIQdJxy1WTPMT61rYfFI7XbjgOcV63oYRnVNcVPHfYYJ\nVxYrUw7busgM6ewA6b7FzZY+ONS5HQ33cV+CNHLCwcFOs7fnNEKqrctLJU07VG9r7BjtWR2b\nNcP81LpOLDbXbjcOeH6xvodhVJ8tVnfaNdjo+LPcsVXFqpTDtpkhfa/11u629kbe2klpkEZ7\nhXlrZ2+8vbUbXPxYWT6/5KH2RiBJQMpuskIaaqxf0mjsLK9Z8eSWKy7c094MJAlI2U1WSMsH\nXHeUL61cdvrVQ2EzkCQgZTdZIXUISBKQsgNSFJAkIGUHpCggSUDKDkhRQJKAlB2QooAkASk7\nIEUBSQJSdkCKApIEpOyAFAUkCUjZASkKSBKQsgNSFJAkIGUHpCggSUDKDkhRQJKAlB2QooAk\nASk7IEUBSQJSdkCKApIEpOyAFAUkCUjZASkKSBKQsgNSFJAkIGUHpCggSUDKDkhRQJKAlB2Q\nooAkASk7IEUBSUqD1GGIEJDsAUkCkgtI9oAkAckFJHtAkoDkApI9IElAcgHJHpAkILmAZA9I\nEpBcQLIHJAlILiDZA5IEJBeQ7AFJApILSPaAJAHJBSR7QJKA5AKSPSBJQHIByR6QJCC5gGQP\nSBKQXECyByQJSC4g2QOSBCQXkOwBSQKSC0j2gCQByQUke0CSgOQCkj0gSUByAckekCQguYBk\nD0gSkFxAsgckCUguINkDkgQkF5DsAUkCkgtI9oAkAckFJHtAkoDkApI9IElAcgHJHpAkILmA\nZA9IEpBcQLIHJAlILiDZA5IEJBeQ7AFJApILSPaAJP2iII28o8cEUrj6yIBkD0gSkFxAsgck\nCUguINkDkgQkF5DsAUkCkgtI9oAkAckFJHtAkoDkApI9IElAcgHJHpAkILmAZA9IEpBcQLIH\nJAlILiDZA5IEJBeQ7AFJApILSPaAJAHJBSR7QJKA5AKSPSBJQHIByR6QJCC5gGQPSBKQXECy\nByQJSC4g2QOSBCQXkOwBSQKSC0j2gCQByQUke0CSgOQCkj0gSUByAckekCQguYBkD0gSkFxA\nsgckCUguINkDkgQkF5DsAUkCkgtI9oAkAckFJHtAkoDkApI9IElAcgHJHpAkILmAZA9IEpBc\nQLIHJAlILiDZA5IEJBeQ7AFJApILSPaAJAHJBSR7QJKA5AKSPSBJQHIByR6QJCC5gGQPSBKQ\nXECyByQJSC4g2QOSBCQXkOwBSQKSC0j2gCQByQUke0CSgOQCkj0gSUByAckekKTJCCkMEQKS\nPSBJQHIByR6QJCC5gGQPSBKQXECyByQJSC4g2QOSBCQXkOwBSQKSC0j2gCQByQUke0CSgOQC\nkj0gSUByAckekCQguYBkD0gSkFxAsgckCUguINkDkgQkF5DsAUkCkgtI9oAkAckFJHtAkoDk\nApI9IElAcgHJHpAkILmAZA9IEpBcQLIHJAlILiDZA5IEJBeQ7AFJApILSPaAJAHJBSR7QJKA\n5AKSPSBJQHIByR6QJCC5gGRvvEM654EdXWq4D/supPaEdbXn3OuQKtecdmjNS9hxwlGbNcP8\n1LoWFo/UbjcOeF6xrodhVNcUN3XcZ5hwZbEy5bCti/oFia9IEl+RsuMrUhSQJCBlB6QoIElA\nyg5IUUCSgJQdkKKAJAEpOyBFAUkCUnZAippAkMLG6gpIISDpgOTvYSClByQdkPw9DKT0gKQD\nkr+HgZQekHRA8vcwkNIDkg5I/h4GUnpA0gHJ38NASg9IOiD5exhI6QFJByR/DwMpPSDpgOTv\nYSClByQdkPw9DKT0gKQDkr+HgZQekHRA8vcwkNIDkg5I/h4GUnpA0gHJ38NASg9IOiD5exhI\n6QFJByR/DwMpPSDpgOTvYSClByQdkPw9DKT0gKQDkr+HgZQekHRA8vcwkNIDkg5I/h4GUnpA\n0gHJ38NASg9IOiD5exhI6QFJByR/DwMpPSDpgOTvYSClByQdkPw9DKT0gKQDkr+HgZQekHRA\n8vcwkNIDkg5I/h4GUnpA0gHJ38NASg9IOiD5exhI6QFJByR/DwMpPSDpgOTvYSClByQdkPw9\nDKT0gKQDkr+HgZQekHRA8vcwkNIDkg5I/h4GUnpA0gHJ38NASg9IOiD5exhI6QFJByR/DwMp\nPSDpgOTvYSClByQdkPw9DKT0gKQDkr+HgZQekHRA8vcwkNIDkg5I/h4GUnpA0gHJ38NASg9I\nOiBV89PWroAUApIOSNX8tLUrIIWApANSNT9t7QpIISDpgFTNT1u7AlIISDogVfPT1q6AFAKS\nDkjV/LS1KyCFgKQDUjU/be0KSCEg6YBUzU9buwJSCEg6IFXz09augBQCkg5I1fy0tSsghYCk\nA1I1P23tCkghIOmAVM1PW7sCUghIOiBV89PWroAUApIOSNX8tLUrIIWApANSNT9t7QpIISDp\ngFTNT1u7AlIISDogVfPT1q6AFAKSDkjV/LS1KyCFgKQDUjU/be0KSCEg6YBUzU9buwJSCEg6\nIFXz09augBQCkg5I1fy0tSsghYCkA1I1P23tCkghIOmAVM1PW7sCUghIOiBV89PWroAUApIO\nSNX8tLUrIIWApANSNT9t7QpIISDpgFTNT1u7AlIISDogVfPT1q6AFAKSDkjV/LS1KyCFgKQD\nUjU/be0KSCEg6YBUzU9buwJSCEg6IFXz09augBQCkg5I1fy0tSsghYCkA1I1P23tCkghIOmA\nVM1PW7sCUghIOiBV89PWroAUApIOSNX8tLUrIIWApANSNT9t7QpIISDpgFTNT1u7AlIISDog\nVfPT1q6AFAKSDkjV/LS1KyCFgKQDUjU/be0KSCEg6YBUzU9buwJSCEg6IFXz09augBQCkg5I\n1fy0tSsghYCkA1I1P23tCkghIOmAVM1PW7sCUghIOiBV89PWroAUApIOSNX8tLUrIIWApANS\nNT9t7QpIISDpgFTNT1u7AlIISDogVfPT1q6AFAKSDkjV/LS1KyCFgKQDUjU/be0KSCEg6YBU\nzU9buwJSCEg6IFXz09augBQCkg5I1fy0tSsghYCkA1I1P23tCkghIOmAVM1PW7sCUghIOiBV\n89PWroAUApIOSNX8tLUrIIWApANSNT9t7QpIISDpgFTNT1u7AlIISDogVfPT1q6AFAKSDkjV\n/LS1KyCFgKQDUjU/be0KSCEg6YBUzU9buwJSCEg6IFXz09augBQCkg5I1fy0tSsghYCkA1I1\nP23tCkghIOmAVM1PW7sCUghIOiBV89PWroAUApIOSNX8tLUrIIWApANSNT9t7QpIISDpgFTN\nT1u7AlIISDogVfPT1q6AFAKSrhdIWy5aUlkBSQJSdpMd0sZlK4GkAlJ2kx3Shu33AUkFpOwm\nO6SyBJIOSNkBqQ3pvsXNlj441LkdDfdxn4fkp61d7X1Ilddr2qH6NWzs6PICd2/WDPNT6zqx\n2Fy73Tjg+cX6HoZRfbZY3WnXYGMw+3SrilUph21b9AuENLSpOE3tmj0lHqiYN7y4sbi2y0kP\nOSDlR9jumOLpbrsHR3uJf3Nq0mUeLRamTtQlINkbx5CkhLd25aPFGWrX3CnRclexYHjxl8UX\nupx01hsSZgx9oHih2+6dL43y/DlTky6zvRhInahLvLWzN47f2klAkoCU3WSHNNRYv6TR2Nle\nAkkCUnaTHdLyAdcd7SWQJCBlN9khjQhIEpCyA1IUkCQgZQekKCBJQMoOSFFAkoCUHZCigCQB\nKTsgRQFJAlJ2QIoCkgSk7IAUBSQJSNkBKQpIEpCyA1IUkCQgZQekKCBJQMoOSFFAkoCUHZCi\ngCQBKTsgRQFJAlJ2QIoCkgSk7IAUBSQJSNkBKQpIEpCyA1IUkCQgZQekKCBJQMoOSFFAkoCU\nHZCigCQBKTsgRQFJAlJ2QIoCkgSk7IAUBSQJSNkBKQpIEpCyA1IUkCQgZQekKCBJQMoOSFFA\nkoCUHZCigCQBKTsgRX181Zc7t/qGL7mHP5v5u2rXEQdFy1tmzhte/OeZf9DlpIf/WpedqgUz\nV3fbffMXR3n+nLcmXeaGmcemTtSlg2erTV+64Sbz6X7j7b0Mozpu5hdqt9/4JdPpTpn5mV6m\nGdnHZ36y066bb7g5+3QXzLwg5bBb+gbpn/5nlz59zlfdw1dXXK12XbYiWt6+4vLhxY0rPt/l\npH9yUbdLjuyKFX/Tbfdtt43y/Mv+S9Jlvr7iM6kTdenCS9WmvznnYvPpLs16pUbtMyu+Xrv9\nG7ebTnfdii/3Ms3IPr/ixk67/vycv8g+3Q0rbkg67u/7Balrl8z5yV48+6Ro55zzxnqEcd9t\nc+78BVwFSPt0QOo9IBGQ+hCQCEh9aPxDIpo0AYmoDwGJqA8BiagP9RXSlouWqG0vfu7Mj1y5\nrXxgQFrbz8tNxLq8hGV518dO/sT3f/Ezjbe6vYZ/5O7CU/fCRfsJaeOylfpHcDmJkJAAAAXn\nSURBVNUlTz573QW7dzWa/ejUbn80j7q/hOW3l23a9s1zR/vjgZO+rq/h2Xc278PBvXDVfkLa\nsP0++REMXbvslE893trWWPxE8/8OTrpfFpf9VR+vNiHr+hKeu2EsRxs3dX0NT9m0l67a318j\ntX4EF137wstf+ejLsuXeD+1pfvzEN9znG5e/0terTcg6v4Q7Bjb80SkX/XgshxsndX4Ndw1c\n/8lzrt6yF665FyA9PjBUlns+vFG2rDvLfbx0TfPD7vPu7uvFJmadX8LNA3/yzAtrPmz/exWT\nps6v4U/P+PPNm68442f9v+ZegLSx9fsKt92zZMmSh9ed7bYLpI1nvdrXi03MOr+Emweab5Bf\nXfrtMR5wHNT1NizLn5+6vv/X3AuQvjfQ+nL60lNPPfXv32t9Tb2t+eHKNX291gSt80vYGHis\n+XjBbWM53fio623Y7A/3wq/V9wKkpwceaX7c2toyuLj5s//8kofK8mf+dxyoa51fwt3L7izL\nl39/41hONz7q/Bo+9fnmL9N3nvqd/l+zn5CGGuuXNBo7y0v/ePur3zrF/x7jNSue3HLFhc3/\nP7h/oP6vKlOlri/hbaf/sHH9sp1jO+G+X7fX8IWPrNy65eqz8//u+aj1E9JyeVN6Rzn0p6f9\n/h8/5De+tHLZ6Vc3f9lX/sNifs9u1Lq+hLtvPePkT/GduNHq+ho+cdlpH71qb/ylBP6IEFEf\nAhJRHwISUR8CElEfAhJRHwISUR8CElEfAhJRHwLSmPS+/ziRLkNAGqNWXl1d/TD7Z2H0Z8gR\n8WWyz0Hp8WLuA12f/bMw+jP6cQSlx4s5Jrn3XAt+69+O2//AD28rTyiKYk5Z/uNv7z9t9i3N\nnfMX3Pm2eWW5/ujXv+VU9zcnwo73HrVh7rQ3nv3T9jOGC4f8v4+9/XVv+eCP/RGty2yc+8sH\nXbvrkoNef/wTzSP+eu60/ef8damuSj0FpDHJ3eHHHzz37m23TzmzfHRJsenh8ttTjr5z/XnF\nn5XlcYe/+4a15frX/M5Xb3nnf9ha2THvwCP+ufGVqSf7Zww3fMhRM2/+ztfe8+aXWkfIZd52\n7L8+c3Lx21du+acDfq8sv16cvHbtwmKtuir1FJDGJLnDi+82Pzv+oLJc7n4WZr/L/ftAi/ff\n2dzxt83Pjvj1V8ryX167qrJjfuH+MtLy4unWM4YLhzxffKr5yeNXP9s6onWZ+8vynuL9zeXp\nM8ry6uNeLsvnf+l0dVXqKSCNSXKHT3efnblf65beVnxyZ7Obiu+Xx792V1nuKP6wdWhlx/wZ\n7u953lp8awSk4UN2vekd3279pyiHITX1lI8XFzc/Xhz+Q6BvW6CuSj0FpDFJ7vB3uM/c7ez+\n98PC97fyRap8sLiidWhlx/x3ug1riy+NgFQ55Lu/XrzpQ197pQrJXeb/FNc0P15SPFc+/18P\nO2DKlGK+uir1FJDGpDpI59wnNVo7flRc3jq0sqMF6ZvFrQpSOKR8dcPF/6k44ucdIR095dMb\nH3jwoPnqqtRTQBqTNKTB4sz2TtnxQiH/8M1T2ys75k9z/w7T6mLdCEiVQ6T/Xny5E6THinOb\nn7zyy/PVVamngDQmxZA+VjTfix35hueay1svfcXveM+BzV/Q/Lj5Bm94x/zmr47K8qTXDbWe\nMVw45AenuX8Y4/HiutYRGtLDxZWl+xbSUeqq1FNAGpNiSJcXV95e/uPUw2/9+8umntXesXa/\n9/3VmkPevLWyY/7Bh9x498XFMv+M4cIhW/c//Ja7v/7+Ax5vHaEh7Tr4rXd896Jjjtn/Oz8b\ncVXqKSCNSTGkZ2ZPba7v+cD+Uw+59pX2jvKuo6a/+eRHy8qO+e/+wdHT33jui+1nDBcO+d8n\nv3nqQSf/mz+i5tdIm+ZNf8sfPH/nr75x84irUk8BaRw1nz+Dus8GpHEUkPbdgDSOAtK+G5DG\nURGkvytCN47ZRNQOSOO1Fx8MDY31LAQkon4EJKI+BCSiPgQkoj4EJKI+BCSiPgQkoj70/wFB\n5hB3yY+uMAAAAABJRU5ErkJggg==",
 510       "text/plain": [
 511        "plot without title"
 512       ]
 513      },
 514      "metadata": {
 515       "image/png": {
 516        "height": 420,
 517        "width": 420
 518       }
 519      },
 520      "output_type": "display_data"
 521     }
 522    ],
 523    "source": [
 524     "ggplot(df, aes(x=intercept_estimate)) + geom_histogram(bins=100) +\n",
 525     "    geom_vline(xintercept=dh$xbr) + \n",
 526     "    scale_x_log10() + scale_y_log10() +\n",
 527     "    theme_bw()\n",
 528     "ggsave(paste0('plot_dhist_', unique(df$op), '.png'))"
 529    ]
 530   }
 531  ],
 532  "metadata": {
 533   "kernelspec": {
 534    "display_name": "R",
 535    "language": "R",
 536    "name": "ir"
 537   },
 538   "language_info": {
 539    "codemirror_mode": "r",
 540    "file_extension": ".r",
 541    "mimetype": "text/x-r-source",
 542    "name": "R",
 543    "pygments_lexer": "r",
 544    "version": "4.0.4"
 545   }
 546  },
 547  "nbformat": 4,
 548  "nbformat_minor": 4
 549 }