30:T5a1,In general, curves closest to the red dot (hypothetical experiment) are most probable, but there may not be *any* curves that pass through the red dot if it represents an outcome that is unlikely for any value of ${unknownParam.name}. Another way to represent the outcome of this hypothetical experiment is shown above (right): a plot showing how likely each value of the unknown ${unknownParam.name} is before (***prior***) and after (***posterior***) the measurement is performed. How do we calculate the posterior probability distribution? It is simply the prior times the likelihood, normalized over parameters, as given by [Bayes theoreom](https://en.wikipedia.org/wiki/Posterior_probability): ${tex.block` {\color{magenta}P(${unknownParam.tex}\mid y_\text{obs}, t_\text{obs})} = \dfrac{ {\color{green}P(${unknownParam.tex})}\, {\color{blue}P(y_\text{obs}\mid ${unknownParam.tex},t)} } { \int d${unknownParam.tex}'\, {\color{green}P(${unknownParam.tex}')}\, {\color{blue}P(y_\text{obs}\mid ${unknownParam.tex}',t)} } `} Go back and try different choices of unknown parameter and prior, and see how they affect these last two plots. Notice how the ***posterior*** can have several peaks, also known as being *multi-modal*, when the frequency is the unknown parameter.31:T69d,## Measuring Information Gain How much information would we gain from our hypothetical experiment? [Information theory](https://en.wikipedia.org/wiki/Information_theory) provides a method of calculating a number for any probability distribution ${tex`p(\theta)`}, called the ***entropy*** ${tex`H(p(\theta))`}, that measures how informative it is: ${tex.block` H(p(\theta)) = -\int p(\theta) \log_2(p(\theta)) \; . `} A flat distribution, where all allowed values are equally likely, has the largest possible entropy, and a narrow peaked distribution, where only a small range of values are likely, has lower entropy. We then define the ***information gain*** of our hypothetical experiment as the difference between the prior entropy and the posterior entropy: see the values shown at the top of the previous plot and watch how they vary with the observation time and measurement result sliders. Notice how the information "gain" can actually be negative in some cases, although this is unlikely to occur in practice. A good experiment provides a large information gain, but when we design our experiment, we only control the observation time ${tex`t_\text{obs}`}, but not the measurement result ${tex`y_\text{obs}`} which is inherently uncertain. We can visualize this situation with the plots below. On the left, a plot of the information gained ${tex`IG(y_\text{obs},t_\text{obs} = ${experiment.tobs.toFixed(1)})`} for all possible measurement outcomes ${tex`y_\text{obs}`}. On the right, a plot showing how likely each measurement outcome is, ${tex`P(y_\text{obs}\mid t_\text{obs} = ${experiment.tobs.toFixed(1)})`}.32:T642,Since the measurement result is uncertain, we calculate an ***expected information gain*** (EIG) as an average of the information gained (IG) with different possible measurement results. However, we need to account for the fact that some measurement results are more likely than others, as shown on the right-hand plot above. If we knew the value of the ${unknownParam.name} ${tex`${unknownParam.tex}`}, the ***likelihood*** would tell us exactly how likely each measurement result is for our hypothetical experiment design with ${tex`t_\text{obs} = ${experiment.tobs.toFixed(1)}`}. However, the whole point of this exercise is to learn about the value of ${tex`${unknownParam.tex}`}, so we cannot do that! Instead, we consider all possible values of ${tex`${unknownParam.tex}`} in the likelihood, weighted by the ***prior***. Specifically, ${tex.block` EIG(t_\text{obs}) = \int dy_\text{obs}\, IG(y_\text{obs}, t_\text{obs}) \, P(y_\text{obs} \mid t_\text{obs}) `} with ${tex.block` P(y_\text{obs} \mid t_\text{obs}) = \int d${unknownParam.tex}\, P(${unknownParam.tex}) P(y_\text{obs}\mid ${unknownParam.tex}, t_\text{obs}) `} where the last integrand is just the prior times the likelihood (similar to the posterior, but without the normalization over ${tex`${unknownParam.tex}`}). Note that while IG can be negative, the weighted average EIG value must always be non-negative. The darkness of the dots on the left-hand plot above shows the calculated weights ${tex`P(y_\text{obs} \mid t_\text{obs}=${experiment.tobs.toFixed(2)})`}.33:T421,designPlot = { const {ny, nxi, y_grid, xi_grid} = grids; const {iy, ixi} = experimentIndices; const IG = Float64Array.from({length: ny}, (_,iy) => S0.IG[iy * nxi + ixi]); const IGmax = d3.max(IG); const EIG = S0.EIG[ixi]; const fmt = d3.format(".2f"); const L = Float64Array.from({length: ny}, (_,iy) => S0.y_marginal[iy * nxi + ixi]); const Lmax = d3.max(L); return Plot.plot({ x: { label: "Measurement result" }, y: { label: "Information gain in bits", domain: [0,1.1*IGmax] }, grid: true, marks: [ Plot.line(IG, { x: (_,iy) => y_grid[iy], y: y => y, stroke: "lightgray" }), Plot.dot(IG, { x: (_,iy) => y_grid[iy], y: y => y, fill: "black", fillOpacity: (_,iy) => L[iy]/Lmax }), Plot.ruleX([y_grid[iy]], { x: x=>x, strokeDasharray: [3,2] }), Plot.text([`EIG = ${EIG.toFixed(2)} bits for measurement time ${xi_grid[ixi].toFixed(1)}`], { text: d=>d, x: y_grid[0], dx: 20, y: 1.1*IGmax, dy: 20, fontSize: 20, textAnchor: "start" }), Plot.frame(), ] }); }34:T625,priorpostPlot = { const { ntheta, theta_grid } = grids; const Pmax = Math.max(d3.max(prior), d3.max(experimentPosterior)); const theta_min = theta_grid[0], theta_max = theta_grid[ntheta-1]; const xtext = frac => theta_min * (1 - frac) + theta_max * frac; const F = d3.format(".1f"); const marks = [ Plot.areaY(prior, { x: (_,i) => theta_grid[i], y: y=>y, fill: "green", fillOpacity:0.25 }), Plot.line(prior, { x: (_,i) => theta_grid[i], y: y=>y, stroke: "green", strokeWidth:2 }), Plot.areaY(experimentPosterior, { x: (_,i) => theta_grid[i], y: y=>y, fill: "magenta", fillOpacity:0.25 }), Plot.line(experimentPosterior, { x: (_,i) => theta_grid[i], y: y=>y, stroke: "magenta", strokeWidth:2 }), Plot.frame(), ]; const postInfoBits = infoBits(experimentPosterior); const legends = [ { text: `Prior entropy ${F(S0.I0)} bits`, color: "green", dx: 8, frac: 0, anchor: "start" }, { text: `Information gain ${F(S0.I0 - postInfoBits)} bits`, color: "black", dx: 0, frac: 0.5, anchor: "middle" }, { text: `Posterior entropy ${F(postInfoBits)} bits`, color: "magenta", dx: -8, frac: 1, anchor: "end" }, ]; legends.forEach((legend,i)=> marks.push( Plot.text([legend], { text: "text", x: d => xtext(d.frac), y: 1.1 * Pmax, dx: legend.dx, dy: 16, fontWeight: "bold", fontSize: 14, fill: "color", textAnchor: legend.anchor }) ) ); return Plot.plot({ x: { label: unknownParam.label }, y: { label: "Prior probability", domain: [0, 1.1 * Pmax] }, grid: true, marks, }); }35:T718,function calculate(prior) { const {log2} = Math; const { ny, ntheta, nxi, decode } = grids; // Calculate the prior information in bits. const I0 = infoBits(prior); // Tabulate the posterior P(theta|y,xi) on the grid. const encode = (iy, ixi) => iy * nxi + ixi; const norm = new Float64Array(ny * nxi); const posterior = Float64Array.from({length: likelihood.length}, (_,n) => { const { itheta, iy, ixi } = decode(n); const P = prior[itheta]; if(Number.isNaN(Number(P))) throw new Error(`n=${n} itheta=${itheta} P=${P}`); const post = likelihood[n] * prior[itheta]; norm[encode(iy, ixi)] += post; return post; }); posterior.forEach((post,n) => { const { iy, ixi } = decode(n); posterior[n] = post / norm[encode(iy, ixi)]; }); // Tabulate the information gain in bits. const IG = new Float64Array(ny * nxi).fill(I0); posterior.forEach((post,n) => { if(post > 0) { const { iy, ixi } = decode(n); IG[encode(iy, ixi)] += post * log2(post); } }); // Tabulate the marginal probability P(y|xi) by integrating P(y|theta,xi) P(theta) over theta. const y_marginal = new Float64Array(ny * nxi); likelihood.forEach((like,n) => { const { iy, itheta, ixi } = decode(n); y_marginal[encode(iy, ixi)] += like * prior[itheta]; }); // Tabulate the weighted iformation gain P(y|xi) IG(y,xi). const wIG = new Float64Array(ny * nxi); IG.forEach((ig,n) => { wIG[n] = y_marginal[n] * ig; }); // Tabulate the expected information gain in bits. const EIG = new Float64Array(nxi); y_marginal.forEach((ym,n) => { const ixi = n % nxi; EIG[ixi] += ym * IG[n]; }); // Find the best design with the maximum EIG. const ixiBest = d3.maxIndex(EIG); return {prior, I0, posterior, IG, y_marginal, wIG, EIG, ixiBest}; }36:T40b,grids = { const {floor, abs} = Math; const ny = 100, ntheta = 64, nxi = 101, nyxi = ny * nxi; const xi_grid = linspace(0, 10, nxi); // theta grid depends on the unknown param const [theta_min, theta_max] = unknownParam.range; const theta_grid = linspace(theta_min, theta_max, ntheta); // y grid depends on sigma and the unknown param let y_lim = (unknownParam.name == "amplitude") ? theta_max : fixedParams.amplitude; y_lim += 3 * sigma; const y_grid = linspace(-y_lim, y_lim, ny); // Function to encode/decode unrolled index n[itheta,iy,ixi] // with n = itheta * ny * nxi + iy * nxi + ixi const encode = ({itheta, iy, ixi}) => itheta * ny * nxi + iy * nxi + ixi; const decode = n => { const itheta = floor(n / nyxi); const iyxi = n % nyxi; const iy = floor(iyxi / nxi); const ixi = iyxi % nxi; return { itheta, iy, ixi }; } const getindex = (grid, x0) => d3.minIndex(grid.map(x => abs(x - x0))); return { ntheta, nxi, ny, theta_grid, xi_grid, y_grid, encode, decode, getindex }; }37:[] 2e:{"title":"Bayesian Optimal Experiment Design","cells":[{"id":589,"value":"# Bayesian Optimal Experiment Design\n\nThis is an interactive introduction to the problem of designing an experiment that is ***optimal***, in the sense of giving us the maximum possible expected gain of information, within the framework of Bayesian statistics. For other introductions, see [wikipedia](https://en.wikipedia.org/wiki/Bayesian_experimental_design) or this [2023 review article](https://arxiv.org/abs/2302.14545).","pinned":false,"hidden":false,"mode":"md"},{"id":0,"value":"## Introduction\n\nImagine that we want to learn as much as possible about how an object is moving but can only make a single observation, or ***experiment***. For concreteness, assume that we know the object has a periodic motion governed by these three unknown parameters. Play with the sliders to get a feel for the full range of possible motions we are considering.","pinned":false,"hidden":false,"mode":"md"},{"id":225,"value":"viewof demoParams = Inputs.form(\n Object.fromEntries(\n paramsConfig.entries().map(\n ([key,value]) => [\n key,\n Inputs.range(value.range, { label:value.label, step:0.01 }),\n ]\n )\n ),\n {\n template: (inputs) => html`\n

${Object.values(inputs)}`\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":248,"value":"modelPlot = {\n const {xi_grid} = grids;\n const {amplitude,omega,offset} = demoParams;\n const [amin,amax] = paramsConfig.get(\"amplitude\").range;\n return Plot.plot({\n x: { label: \"Observation time t\" },\n y: { label: \"Displacement y(t)\", domain: [-amax,+amax] },\n grid: true,\n marks: [\n Plot.line(xi_grid, { x: x=>x, y: x=>y_mean(x, demoParams), stroke: \"blue\" }),\n Plot.frame(),\n ]\n })\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":270,"value":"The corresponding mathematical model is:\n${tex.block`\ny(t) = {\\color{red}A} \\sin\\left[ {\\color{red}\\omega} (t - {\\color{red}t_0})\\right]\n`}\nwhere the three unknown parameters are highlighted in red.\n\nTo simplify the calculations and the interpretation of our results, let's assume that we already know the values of two of these parameters, so we only have one unknown parameter to determine with our experiment. Pick which parameter is unknown below:","pinned":false,"hidden":false,"mode":"md"},{"id":283,"value":"viewof unknownParam = {\n const options = [ ...paramsConfig.entries() ].map(([key,value]) => [value.label, value]);\n return Inputs.radio(new Map(options), {label: \"Unknown parameter\", value: options[1][1]});\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":353,"value":"The two fixed parameters are\n${[...paramsConfig.values()].filter(({name}) => name != unknownParam.name).map(\n ({name,label}) => `${label} = ${fixedParams[name]}`).join(\" and \")},\n which you can adjust with these sliders:","pinned":false,"hidden":false,"mode":"md"},{"id":507,"value":"viewof fixedParams = Inputs.form(\n Object.fromEntries(\n [...paramsConfig.values()]\n .filter(\n (param) => param.name != unknownParam.name\n )\n .map(\n (param) => [param.name, Inputs.range(param.range, { label:param.label, step:0.01 })]\n )\n ),\n {\n template: (inputs) => html`\n

${Object.values(inputs)}`\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":351,"value":"The goal of our experiment is now to learn which of the curves below best describes the motion of the object:","pinned":false,"hidden":false,"mode":"md"},{"id":615,"value":"priorCurves = modelsGridPlot(prior)","pinned":false,"hidden":false,"mode":"ojs"},{"id":383,"value":"## Prior Assumptions\n\nBefore we can design an optimal experiment to determine the unknown ${unknownParam.name} ${tex`${unknownParam.tex}`}, we must specify its range of possible values and how likely different values within this range are *before we make any measurement*. These constitute our ***prior*** ${tex`P(${unknownParam.tex})`} on the ${unknownParam.name} ${tex`${unknownParam.tex}`}. Choose one of the following priors to proceed (and note how your selection is represented in the previous and following plot):","pinned":false,"hidden":false,"mode":"md"},{"id":393,"value":"viewof priorSel = Inputs.radio(priorsConfig, {label: \"Prior\", value: [...priorsConfig.values()][1] })","pinned":false,"hidden":false,"mode":"ojs"},{"id":404,"value":"priorPlot = {\n const { theta_grid } = grids;\n const Pmax = d3.max(prior);\n const dtheta = theta_grid[1] - theta_grid[0];\n return Plot.plot({\n x: { label: unknownParam.label },\n y: { label: \"Prior probability\", domain: [0, 1.1 * Pmax] },\n grid: true,\n marks: [\n Plot.areaY(prior, { x: (_,i) => theta_grid[i], y: y=>y, fill: \"green\", fillOpacity:0.25 }),\n Plot.line(prior, { x: (_,i) => theta_grid[i], y: y=>y, stroke: \"green\", strokeWidth:2 }),\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":435,"value":"## Measurement Uncertainty and Likelihood\n\nWe need one final ingredient before we can design an optimal experiment: measurement uncertainty. Without uncertainty, almost any measurement would reveal the motion perfectly. Experiment design is all about optimizing what we can learn in the presence of uncertainty, given our prior beliefs.\n\nWe quantify measurement uncertainty as the ***likelihood*** ${tex`P(y_\\text{obs}\\mid ${unknownParam.tex},t)`} of measuring a value ${tex`y_\\text{obs}`} given an assumed value of the unknown ${unknownParam.name} ${tex`${unknownParam.tex}`} and our chosen measurement time ${tex`t_\\text{obs}`}. For this example, we assume that measurement errors\n${tex.block`\n\\Delta y \\equiv y_\\text{obs} - y(t_\\text{obs})\n`}\nhave a bell-shaped or Gaussian distribution with a known spread that you can set below:","pinned":false,"hidden":false,"mode":"md"},{"id":526,"value":"viewof sigma = Inputs.radio(\n new Map([[\"small\", 0.2],[\"medium\", 0.4],[\"large\", 0.8]]),\n { label: \"Measurement uncertainty\", value: 0.4 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":504,"value":"Here we visualize the likelihood for ${unknownParam.name} ${tex`${unknownParam.tex} = ${likelihoodTheta.toFixed(2)}`} (adjustable with the slider below) where brighter regions show the most likely measured values ${tex`y_\\text{obs}`} at each measurement time ${tex`t`}. Notice how larger uncertainties blur the underlying model vertically.","pinned":false,"hidden":false,"mode":"md"},{"id":524,"value":"viewof likelihoodTheta = {\n const { theta_grid, ntheta } = grids;\n return Inputs.range([theta_grid[0], theta_grid[ntheta-1]], {label: \"Likelihood for \" + unknownParam.label, format: x => x.toFixed(2) });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":456,"value":"likelihoodPlot = {\n const { ny, nxi, theta_grid, getindex } = grids;\n const itheta = getindex(theta_grid, likelihoodTheta);\n const nyxi = ny * nxi, n0 = itheta * nyxi;\n return rasterPlot(likelihood.slice(n0, n0 + nyxi), \"Likelihood\");\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":596,"value":"## Outcome of a Hypothetical Experiment\n\nWe now have all of the ingredients needed to analyze the outcome of a hypothetical experiment, which consists of a chosen observation time ${tex`${grids.xi_grid[0]} \\le t_\\text{obs} \\le ${grids.xi_grid[grids.nxi-1]}`} and a resulting (uncertain) measurement result ${tex`y_\\text{obs}`}:","pinned":false,"hidden":false,"mode":"md"},{"id":601,"value":"viewof experiment = {\n const { y_grid, xi_grid, ny, nxi } = grids;\n return Inputs.form({\n tobs: Inputs.range([xi_grid[0], xi_grid[nxi-1]], {label: \"Observation time\", format: x => x.toFixed(1), step: 0.1 }),\n yobs: Inputs.range([y_grid[0], y_grid[ny-1]], {label: \"Measurement result\", format: x => x.toFixed(1), step: 0.1}),\n },\n {\n template: (inputs) => html`\n

${Object.values(inputs)}`\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":610,"value":"Our hypothetical measurement would narrow down the possible curves describing the unknown motion, as represented below (left) with darker curves being more probable. Compared with the similar plot above of possible curves *before* any measurement, which represents our ***prior*** knowledge, this new plot represents our ***posterior*** knowledge.","pinned":false,"hidden":false,"mode":"md"},{"id":913,"value":"priorpostPlots = html`\n

\n ${priorpostCurves} ${priorpostPlot}\n

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":688,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":706,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":939,"value":"igPlots = html`\n

\n ${designPlot} ${ymargPlot}\n

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":783,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":892,"value":"Notice how the most probable measurement results (peaks in this plot) generally correspond to the smallest information gains (valleys in the previous plot). This is a consequence of how we define ***information***: a more certain outcome is already expected so provides less useful information.","pinned":false,"hidden":false,"mode":"md"},{"id":822,"value":"## What is the Optimal Experiment?\n\nWe are finally ready to reveal the optimal experiment design, which consists of picking at what time ${tex`t_\\text{obs}`} to perform a measurement in order to achieve the maximum possible ***expected information gain***. We can read this straight off the plot below of information gain versus measurement time, finding that the best value is ${tex`t_\\text{obs} = ${grids.xi_grid[S0.ixiBest].toFixed(2)}`} where the expected information gain is ${S0.EIG[S0.ixiBest].toFixed(2)} bits (indicated with a blue dot below).","pinned":false,"hidden":false,"mode":"md"},{"id":824,"value":"eigPlot = {\n const { xi_grid } = grids;\n const { EIG, ixiBest } = S0;\n const xiBest = xi_grid[ixiBest];\n return Plot.plot({\n x: { label: \"Measurement time\" },\n y: { label: \"Expected information gain in bits\", domain: [0,1.1*d3.max(EIG)] },\n grid: true,\n marks: [\n Plot.line(EIG, { x: (_,i) => xi_grid[i], y: y=>y, stroke: \"black\" }),\n Plot.ruleX([xiBest], { x: x=>x, stroke: \"blue\" }),\n Plot.dot([xiBest], { x: x=>x, y: EIG[ixiBest], r: 4, fill: \"blue\" }),\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":861,"value":"## Next Steps\n\nSuppose we perform the optimal experiment but still need a more accurate estimate of the parameter ${tex`${unknownParam.tex}`}, what next? We can simply repeat the entire process outlined above, but replacing our original ***prior*** with the computed ***posterior*** for the actual measurement result we obtain. This recipe can then be repeated as many times as needed, to design a sequence of experiments that optimally zero in on the true value of the unknown parameter ${tex`${unknownParam.tex}`}.","pinned":false,"hidden":false,"mode":"md"},{"id":208,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":785,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":865,"value":"ymargPlot = {\n const {ny, nxi, y_grid, xi_grid} = grids;\n const {iy, ixi} = experimentIndices;\n const L = Float64Array.from({length: ny}, (_,iy) => S0.y_marginal[iy * nxi + ixi]);\n const Lmax = d3.max(L);\n return Plot.plot({\n x: { label: \"Measurement result\" },\n y: { label: \"Relative probability of measured values\", domain: [0,1.1*Lmax] },\n grid: true,\n marks: [\n Plot.line(L, { x: (_,iy) => y_grid[iy], y: y => y, stroke: \"gray\" }),\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":700,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":668,"value":"priorpostCurves = {\n const {tobs,yobs} = experiment;\n const opts = { stroke: \"red\", strokeWidth: 3, strokeOpacity: 0.5, strokeDasharray: [ 6, 2 ], clip:true };\n const extras = [\n Plot.ruleX([tobs], opts),\n Plot.ruleY([yobs], opts),\n Plot.dot([{x:tobs, y:yobs}], { x:\"x\", y:\"y\", fill:\"red\", r:6, clip:true }),\n ];\n return modelsGridPlot(experimentPosterior, { extras });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":769,"value":"","pinned":false,"hidden":false,"mode":"html"},{"id":124,"value":"function rasterPlot(values, label) {\n const { nxi, ny, xi_grid, y_grid } = grids;\n const vmin = d3.min(values);\n const ctype = vmin < 0 ? \"diverging\" : \"sequential\";\n const scheme = vmin < 0 ? \"RdBu\" : \"Magma\";\n return Plot.plot({\n x: {label: \"Measurement time t\"},\n y: {label: \"Measured value y\", domain: [-2.5, 2.5]},\n color: {label: label, legend: true, scheme, type: ctype },\n marks: [\n Plot.rect([0], { x1: xi_grid[0], x2: xi_grid[nxi-1], y1: -2.5, y2: +2.5, fill: \"#000\" }),\n Plot.raster(values, {\n width: nxi, height: ny, x1: xi_grid[0], x2: xi_grid[nxi-1], y1: y_grid[0], y2: y_grid[ny-1],\n clip: true, interpolate: \"none\"\n }),\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":307,"value":"function modelsGridPlot(prob, { extras } = {}) {\n const unknown = unknownParam.name;\n const { xi_grid, theta_grid } = grids;\n const [amin,amax] = paramsConfig.get(\"amplitude\").range;\n const pmax = d3.max(prob);\n const marks = models_grid.map(\n (yvalues, i) => Plot.line(yvalues, {\n x: (_,i) => xi_grid[i], y: y => y, stroke: \"blue\",\n strokeOpacity: prob[i * models_grid_downsample] / pmax\n }));\n marks.push(extras);\n marks.push(Plot.frame());\n return Plot.plot({\n x: { label: \"Observation time t\" },\n y: { label: \"Displacement y(t)\", domain: [-amax,+amax] },\n grid: true,\n marks\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":732,"value":"experimentInfo = {\n const { log2 } = Math;\n return -d3.sum(experimentPosterior, prob => prob * log2(prob));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":659,"value":"experimentPosterior = {\n const { tobs, yobs } = experiment;\n const { ntheta, ny, nxi } = grids;\n const { iy, ixi } = experimentIndices;\n const nyxi = ny * nxi, n0 = iy * nxi + ixi;\n return Float64Array.from({length: ntheta}, (_,itheta) => S0.posterior[n0 + itheta*nyxi]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":726,"value":"experimentIndices = {\n const { tobs, yobs } = experiment;\n const { y_grid, xi_grid, getindex } = grids;\n const iy = getindex(y_grid, yobs), ixi = getindex(xi_grid, tobs);\n return { iy, ixi };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":101,"value":"S0 = calculate(prior)","pinned":false,"hidden":false,"mode":"ojs"},{"id":98,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":739,"value":"function infoBits(prob) {\n const { log2 } = Math;\n return -d3.sum(prob, p => p > 0 ? p * log2(p) : 0);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":25,"value":"likelihood = {\n const { ntheta, nxi, ny, theta_grid, xi_grid, y_grid, encode, decode } = grids;\n const unknown = unknownParam.name;\n const params = Object.assign({}, fixedParams); // shallow copy\n const L = new Float64Array(ny * ntheta * nxi);\n for(let itheta=0; itheta < ntheta; ++itheta) {\n params[unknown] = theta_grid[itheta];\n for(let ixi=0; ixi < nxi; ++ixi) {\n const mu = y_mean(xi_grid[ixi], params);\n let norm = 0.0;\n for(let iy=0; iy < ny; ++iy) {\n const like = Math.exp(-0.5 * ((y_grid[iy] - mu) / sigma) ** 2);\n const n = encode({itheta, iy, ixi});\n L[n] = like;\n norm += like;\n }\n // Second pass to normalize for this value of (theta, xi)\n for(let iy=0; iy < ny; ++iy) {\n const n = encode({itheta, iy, ixi});\n L[n] /= norm;\n }\n }\n }\n return L;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":396,"value":"prior = {\n const { theta_grid, ntheta } = grids;\n const lo = theta_grid[0], hi = theta_grid[ntheta - 1];\n const unnormalized = theta_grid.map(theta => priorSel.P(theta,lo,hi));\n const norm = d3.sum(unnormalized);\n return unnormalized.map(P => P / norm);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":391,"value":"priorsConfig = new Map([\n [\"flat\", { P: (theta) => 1, label: \"All values equally likely\" } ],\n [\"central\", { P: (theta,lo,hi) => 1 + Math.cos((theta - (lo+hi)/2) * Math.PI/((hi-lo)/2)), label: \"Central values more likely\" }],\n [\"low\", { P: (theta,lo,hi) => 1 - (theta - lo) / (hi - lo), label: \"Low values more likely\" }],\n [\"high\", { P: (theta,lo,hi) => (theta - lo) / (hi - lo), label: \"High values more likely\" }]\n])","pinned":false,"hidden":false,"mode":"ojs"},{"id":309,"value":"models_grid = {\n const {ntheta, nxi, xi_grid, theta_grid} = grids;\n const unknown = unknownParam.name;\n const params = Object.assign({}, fixedParams); // shallow copy\n return Array.from(theta_grid).filter((_,i) => i % models_grid_downsample == 0).map((theta) => {\n params[unknown] = theta;\n return Float32Array.from({length: nxi}, (_,i) => y_mean(xi_grid[i], params));\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":427,"value":"models_grid_downsample = 2","pinned":false,"hidden":false,"mode":"ojs"},{"id":232,"value":"paramsConfig = new Map([\n [ \"amplitude\", { range: [0.75, 1.25], label: \"Amplitude A\", name: \"amplitude\", tex:\"A\" } ],\n [ \"omega\", { range: [0.5,2], label: \"Frequency ω\", name: \"omega\", tex: \"\\\\omega\" } ],\n [ \"offset\", { range: [-0.5,0.5], label: \"Offset t0\", name: \"offset\", tex: \"t_0\" } ]\n])","pinned":false,"hidden":false,"mode":"ojs"},{"id":473,"value":"//test_decode()","pinned":true,"hidden":false,"mode":"ojs"},{"id":57,"value":"function test_decode() {\n const { ntheta, nxi, ny, encode, decode } = grids;\n const eq = d3.cross(d3.range(ntheta), d3.range(ny), d3.range(nxi)).map(([i,j,k],n) => {\n const { itheta, iy, ixi } = decode(n);\n const n2 = encode({itheta, iy, ixi});\n return (n == n2) && (i == itheta) && (j == iy) && (k == ixi);\n });\n return eq.every(x => x);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":48,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"function y_mean(t, {amplitude=1, omega=1, offset=0}={}) {\n return amplitude * Math.sin(omega * (t - offset));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"import { linspace } from \"@dkirkby/jumpy\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W37","stdlib":"1","titling":"auto"}