2b:null 2c:["$","title",null,{"children":"Approximate Bayesian Computation | eliaskal | Observable"}] 31:[] 2f:{"id":"f8bc463bd88390a1","slug":"approximate-bayesian-computation","title":"Approximate Bayesian Computation","thumbnail":null,"access_level":"public","create_time":"$D2018-02-07T16:24:55.707Z","update_time":"$D2018-02-10T19:57:48.971Z","view_time":null,"first_public_time":"$D2018-02-08T17:14:09.779Z","fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"b131ffdd9aa4cfd1","login":"eliaskal","name":"","avatar_url":"https://avatars-next.observableusercontent.com/78da7d75adb465c46ef217444817b410dca02778e88a0e79d8071ef3b7bf7c2b","tier":"starter_2024"},"version":226,"publish_version":226,"role":null,"authored":false,"authors":[{"type":"individual","id":"41a188d61429a4e2","login":"eliaskal","name":"","avatar_url":"https://avatars-next.observableusercontent.com/78da7d75adb465c46ef217444817b410dca02778e88a0e79d8071ef3b7bf7c2b","tier":null}],"members":[{"type":"individual","id":"41a188d61429a4e2","login":"eliaskal","name":"","avatar_url":"https://avatars-next.observableusercontent.com/78da7d75adb465c46ef217444817b410dca02778e88a0e79d8071ef3b7bf7c2b","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"Approximate Bayesian Computation","cells":[{"id":0,"value":"md`# Approximate Bayesian Computation\n\nThe following are some notes regarding [Approximate Bayesian Computation (ABC)](https://en.wikipedia.org/wiki/Approximate_Bayesian_computation), which is a simple method for inferring the probability distribution of an unknown parameter.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":154,"value":"md`\n## Example Markov model\n\nLet's consider a Markov model, where there are two states, A and B, and the transition probability from either state to the other is *θ* (the example is based on [this article](https://en.wikipedia.org/wiki/Approximate_Bayesian_computation)).\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":162,"value":"dot`\ndigraph {\n\t{\n\t\tnode [shape = circle]\n\t\tA\n\t\tB\n\t}\n\tA -> B [label = \"θ\"];\n\tB -> A [label = \"θ\"];\n\tA -> A [label = \"1 - θ\"];\n\tB -> B [label = \"1 - θ\"];\n}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":160,"value":"md`\nThis model could be used to produce sequences of states of arbitrary size. Let’s create a function simulating such a sequence.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":8,"value":"function toggle(state) {\n return (state == \"A\" ? \"B\" : \"A\");\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":19,"value":"function markov(n, theta) {\n let arr = Array(n);\n // choose first item randomly\n arr[0] = (Math.random() < 0.5) ? \"A\" : \"B\";\n for (let i = 1; i < n; i++) {\n arr[i] = (Math.random() < theta) ? toggle(arr[i-1]) : arr[i-1];\n }\n return arr;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":21,"value":"md`Let's try this function with some examples:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":11,"value":"markov(10, 0.5)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":27,"value":"markov(10, 0.9)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":30,"value":"markov(15, 0.1)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":33,"value":"md`Let's also create a function that counts how many transitions there are in a given sequence.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":35,"value":"function countTrans(arr) {\n let cnt = 0;\n for (let i = 1; i < arr.length; i++) {\n if (arr[i] != arr[i-1])\n cnt++;\n }\n return cnt;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":38,"value":"md`For example:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":40,"value":"countTrans([\"A\", \"B\", \"B\", \"B\", \"A\", \"A\", \"A\", \"B\", \"B\", \"A\"])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":43,"value":"md`## Parameter inference\n\nLet's now suppose that we do not know the *θ* parameter of the Markov model and wish to infer it from a given sequence. Let's suppose that the true *θ* value is ${thetaTrue} and the given sequence, of size ${dataSeq.length}, is the following:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":51,"value":"dataSeq = markov(dataSeqSize, thetaTrue)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":56,"value":"md`\nOne way we can proceed in order to infer the *θ* parameter from the data is with Approximate Bayesian Computation (ABC). We will take random choices of theta, generate sequences and see if the number of transitions is similar to the number of transitions of the given data. If it is, we will keep the value of *θ* chosen in an array and move on to the next choice of *θ*.\n\nSo, let's create a function to generate an array of possible *θ* values from which the given number of transitions could have occured. Our initial guesses for the *θ* value will be taken from a uniform distribution over [0,1]. A tolerance argument will also be used, in order to control which *θ* values to allow to pass.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":58,"value":"function abcInf(nSamples, targetSeq, tolerance) {\n let samples = Array(nSamples);\n let c = 0\n while (c < nSamples) {\n let theta = Math.random();\n let testSeq = markov(targetSeq.length, theta);\n if (Math.abs(countTrans(testSeq) - countTrans(targetSeq)) <= tolerance) {\n samples[c] = theta;\n c++;\n }\n }\n return samples;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":63,"value":"md`Let's try this function for producing some *θ* samples.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":65,"value":"abcInf(10, dataSeq, 1)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":67,"value":"md`\nWe can see that there are some values close to the true value of ${thetaTrue}, but there are also some very different values, since they could also have produced the target number of transitions.\n\nLet's take several *θ* samples and make a histogram of them, to see the approximate distribution of the *θ* values.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":69,"value":"samples = abcInf(nSamples, dataSeq, 1)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":76,"value":"chart = {\n const svg = d3.select(DOM.svg(width, height));\n \n const bar = svg.append(\"g\")\n .attr(\"fill\", \"steelblue\")\n .selectAll(\"rect\")\n .data(bins)\n .enter().append(\"rect\")\n .attr(\"x\", d => x(d.x0) + 1)\n .attr(\"width\", d => x(d.x1) - x(d.x0) - 1)\n .attr(\"y\", d => y(d.length))\n .attr(\"height\", d => y(0) - y(d.length));\n\n svg.append(\"g\")\n .call(xAxis);\n \n svg.append(\"g\")\n .call(yAxis);\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":145,"value":"md`\nYou can experiment with more values of \\`thetaTrue\\`, \\`dataSeqSize\\` and \\`nSamples\\`, using the sliders below, to see their effect to the histogram.\nThe larger the dataSeqSize, the narrower the distribution, since larger target sequences make us more confident about the value of *θ*. Larger number of samples make the distribution smoother, as it more closely approximates the real distribution.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":125,"value":"viewof thetaTrue = html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":129,"value":"thetaTrue","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":132,"value":"viewof dataSeqSize = html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":138,"value":"dataSeqSize","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":140,"value":"viewof nSamples = html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":142,"value":"nSamples","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":208,"value":"md`\n## Why use a sampling procedure\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":201,"value":"md`\nLet's take the mean value of the produced samples:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":203,"value":"sampleMean = {\n let sum = 0;\n for (let i = 0; i < samples.length; i++) {\n sum += samples[i];\n }\n return sum / samples.length;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":205,"value":"md`\nThe mean value is close to the true value of ${thetaTrue}. As a comparison, let's take the number of transitions in our given sequence and divide it with the total number of possible transitions, which is simply the length of the sequence minus one:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":207,"value":"naiveMean = countTrans(dataSeq) / (dataSeq.length - 1)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":213,"value":"md`\nWe can see that the value we get when we just count the number of transitions in the given sequence is very close to the value we get from the sampling procedure.\n\nSo why bother going through ABC, when we could have reached a similar result in a much simpler way?\n\nThe reason is that ABC approximates the *distribution* of the *θ* value, given an observed sequence, which gives us a lot more information than just taking the mean value. One important piece of information that a distribution gives us is how much concentrated it is around its mean. This is a measure of how confident one can be that the computed mean value is close to the true value of the parameter.\n\nThis amount of confidence is very important in order to make judgements about the computed parameters of a model and whether the available data are sufficient to make these judgements. This is related to the concept of [statistical power](https://www.statisticsdonewrong.com/power.html), in statistics.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":189,"value":"md`\n## Appendix\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":122,"value":"bins = d3.histogram()\n .domain(x.domain())\n .thresholds(x.ticks(40))\n (samples)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":80,"value":"x = d3.scaleLinear()\n .domain(d3.extent([0, 1])).nice()\n .range([margin.left, width - margin.right])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":85,"value":"y = d3.scaleLinear()\n .domain([0, d3.max(bins, d => d.length)]).nice()\n .range([height - margin.bottom, margin.top])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":89,"value":"xAxis = g => g\n .attr(\"transform\", `translate(0,${height - margin.bottom})`)\n .call(d3.axisBottom(x).tickSizeOuter(0))\n .call(g => g.append(\"text\")\n .attr(\"x\", width - margin.right)\n .attr(\"y\", -4)\n .attr(\"fill\", \"#000\")\n .attr(\"font-weight\", \"bold\")\n .attr(\"text-anchor\", \"end\")\n .text(samples.x))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":92,"value":"yAxis = g => g\n .attr(\"transform\", `translate(${margin.left},0)`)\n .call(d3.axisLeft(y))\n .call(g => g.select(\".domain\").remove())","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":87,"value":"height = 300","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":82,"value":"margin = ({top: 20, right: 20, bottom: 30, left: 40})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":74,"value":"d3 = require(\"https://d3js.org/d3.v5.min.js\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":164,"value":"dot = {\n const Viz = await require(\"https://unpkg.com/viz.js@1.8.0/viz-lite.js\").catch(() => window.Viz);\n return function(strings) {\n let string = strings[0] + \"\", i = 0, n = arguments.length;\n while (++i < n) string += arguments[i] + \"\" + strings[i];\n string = Viz(string);\n return html([string.slice(string.indexOf(\"