31:T42f,distForecasts = { return forecasts.map(forecast => { const kernelWidth = d3.deviation(data.filter(d => d.t < forecast.baseT), d => d[vField]) * kw// * 2 //const kernelWidth = d3.deviation(forecast.nexts, d => d[vField]) * kw * 2 const VW = forecast.nexts.map(d => [d[vField], d.w]) const kernels = kdeFunctions(VW, {kernelWidth}) const baseDist = kernels.pdf const area = areaUnderFunc(baseDist, kernels.relDomain) const dist = d => baseDist(d) / area const mean = d3.mean(forecast.nexts, d => d[vField]) const std = d3.deviation(forecast.nexts, d => d[vField]) const normal = gaussian(mean, std**2) const normalDist = d => normal.pdf(d) let distExtent = d3.extent(forecast.nexts, d => d[vField]) distExtent = [distExtent[0]-std*3, distExtent[1]+std*3] const distPeaks = getPeaks(dist, distExtent) return { baseT: forecast.baseT, tp: forecast.tp, actual: forecast.actual, point: forecast[vField], dist, normalDist, distPeaks, other: {nexts: forecast.nexts, distExtent} } }) }32:T432,## Distribution Prediction Scores It is useful to have a benchmark distributional forecasting method to compare to the KDE method (which we use to generate the shaded forecasting area). We have chosen a normal distribution method for this purpose: instead of estimating a distribution using kernels over the neighbors' futures, this method takes the mean and standard deviation of the neighbors' futures and uses these for a normal distribution. A distributional forecast is a (/probabilistic) forecast which is a distribution of values, in contrast to a point forecast which is just a single value. Distributional forecasts need [different techniques](https://arxiv.org/pdf/1709.04743.pdf) to compare their accuracy. A popular technique is the **logarithmic score**, which we have chosen here. A comparison of the forecast scores of the KDE method and the normal distribution method shows whether the principles behind the KDE method are useful in practice. Is a potentially multi-model, nuanced summary of the forecasts useful? Or would the normal distribution suffice?33:T421,scoreResults = distForecasts.map(d => { const minPdf = 0.0000000000001 // Prevent infinity in log // const distPdf = d.dist(d.actual) // const normalPdf = d.normalDist(d.actual) // // const distScore = -Math.log(distPdf == 0 ? minPdf : distPdf) // // const normalScore = -Math.log(normalPdf == 0 ? minPdf : normalPdf) // const distScore = -Math.log(Math.max(distPdf, minPdf)) // const normalScore = -Math.log(Math.max(normalPdf, minPdf)) const scoreFunc = negativeLogScore//naiveNegativeLinearScore const distScore = scoreFunc(d.dist, d.actual) const normalScore = scoreFunc(d.normalDist, d.actual) const pointError = Math.abs(d.point - d.actual) const distPeakErrors = d.distPeaks.map(peak => Math.abs(peak[0]-d.actual)) const distPeakErrorWorst = d3.max(distPeakErrors) const distPeakErrorBest = d3.min(distPeakErrors) const distPeakErrorMean = d3.mean(distPeakErrors) return {baseT: d.baseT, tp: d.tp, distScore, normalScore, pointError, distPeakErrors, distPeakErrorWorst, distPeakErrorBest, distPeakErrorMean} })34:T4b6,{ const mean = d3.mean(plotForecast.other.nexts, d => d[vField]) const std = d3.deviation(plotForecast.other.nexts, d => d[vField]) let plotExtent = d3.extent(distForecasts, forecast => d3.extent(forecast.other.nexts, d => d[vField]))//[mean-4*std, mean+4*std] plotExtent = [d3.min(plotExtent[0])-std*3, d3.max(plotExtent[1])+std*3] // const globalMean = d3.mean(data, d => d[vField]) // const globalStd = d3.deviation(data, d => d[vField]) // const plotExtent = [globalMean - globalStd*6, globalMean + globalStd*6] const normPoints = funcPoints(plotForecast.normalDist, plotExtent, 300) const distPoints = funcPoints(plotForecast.dist, plotExtent, 300) return Plot.plot({ marginLeft: 70, height: 260, width: width, //x: {domain: plotExtent}, marks: [ Plot.ruleY([0]), Plot.ruleX([plotForecast.actual], {stroke: "green", strokeDasharray:"2,2"}), Plot.lineY(normPoints, {x: d => d[0], y: d => d[1], stroke: "grey"}), Plot.lineY(distPoints, {x: d => d[0], y: d => d[1], stroke: "blue"}), Plot.dot(plotForecast.other.nexts, {x: d => d[vField], y: 0}), Plot.dot(plotForecast.distPeaks, {x: d => d[0], y: d =>d[1]}) ] }) }35:T4c6,{ const tpPlot = Plot.plot({ height: 260, width: width/2, y: {label: "MAE"}, marks: [ Plot.lineY(horizonPeakResults, {x: "tp", y: "meanPointError", stroke: "grey", strokeDasharray: "1,3"}), Plot.dot(horizonPeakResults, {x: "tp", y: "meanPointError", stroke: "grey"}), Plot.lineY(horizonPeakResults, {x: "tp", y: "meanDistPeakErrorBest", stroke: "blue", strokeDasharray: "1,3"}), Plot.dot(horizonPeakResults, {x: "tp", y: "meanDistPeakErrorBest", stroke: "blue"}), Plot.lineY(horizonPeakResults, {x: "tp", y: "meanDistPeakErrorWorst", stroke: "red", strokeDasharray: "1,3"}), Plot.dot(horizonPeakResults, {x: "tp", y: "meanDistPeakErrorWorst", stroke: "red"}), // Plot.lineY(horizonPeakResults, {x: "tp", y: "meanDistPeakErrorMean", stroke: "green", strokeDasharray: "1,3"}), // Plot.dot(horizonPeakResults, {x: "tp", y: "meanDistPeakErrorMean", stroke: "green"}), ] }) const tpPlot2 = Inputs.table(horizonPeakResults, { maxWidth: width/2-2, width: {tp: 20}, height: 260, format: { diff: d => htl.html`

${d.toFixed(3)}

` }}) return plotFlex([tpPlot, tpPlot2]) }36:[] 2f:{"title":"EpiForecast: Distributional Accuracy","cells":[{"id":0,"value":"# EpiForecast: Distributional Accuracy\n\nAccompaniment to the [EpiForecast](https://episphere.github.io/forecast/) tool. EpiForecast is an exploratory forecast tool which produces interpretable forecasts using interactive visualization and a simple and intuitive forecasting technique (based off empirical dynamic modeling). The resulting forecasts are designed to augment point forecasts with visualization and interactivity, meaning it is difficult to use traditional techniques to compare the results with other forecasting methods. However, for a given dataset it is useful to compare the forecasts with some basic benchmarks, and so this notebook provides two things:\n\n* A comparison of distributional accuracy of our kernel-density estimation (KDE) variant with a typical normal distribution method.\n* A comparison of the point accuracy of the simplex point forecast method and the best and worst peaks from the KDE distribution.\n\nMore information is provided throughout the notebook.","pinned":false,"hidden":false,"mode":"md"},{"id":32,"value":"## Data","pinned":false,"hidden":false,"mode":"md"},{"id":11,"value":"import {dataAllCause as data} from \"@siliconjazz/cdc-weekly-mortality\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":25,"value":"viewof vField = Inputs.select(Object.keys(data[0]), {label: \"Value Field\", value: \"deaths\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":34,"value":"## Forecasting\n\nThe forecasting parameters. For more explanation, see the [sister notebook](https://observablehq.com/@siliconjazz/edm-interpretable-forecasting).","pinned":false,"hidden":false,"mode":"md"},{"id":40,"value":"viewof E = Inputs.range([0, 32], {label: \"E\", step: 1, value: 10})","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"viewof tp = Inputs.range([0, 32], {label: \"Forecast Horizon\", step: 1, value: 32})","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"viewof nn = Inputs.range([0, 32], {label: \"Num. Neighbors\", step: 1, value: 8})","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"forecasts = {\n const dataByT = d3.index(data, d => d.t)\n const forecasts = fc.simplex(data, vField, tp, E, nn, 1)\n forecasts.forEach(forecast => {\n forecast.actual = dataByT.get(forecast.t)?.[vField]\n })\n return forecasts.filter(d => d.actual != null)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":71,"value":"viewof kw = Inputs.range([0, 4], {label: \"Rel. Kernel Width\", value: 1})","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":85,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":496,"value":"Lower is better!","pinned":false,"hidden":false,"mode":"html"},{"id":94,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":467,"value":"### Overall","pinned":false,"hidden":false,"mode":"md"},{"id":477,"value":"{\n const overall = [\n {dist: \"kernel\", meanScore: d3.mean(scoreResults, d => d.distScore)},\n {dist: \"normal\", meanScore: d3.mean(scoreResults, d => d.normalScore)},\n ]\n return Plot.plot({\n width: 400,\n height: 300,\n grid: true,\n marks: [\n Plot.barY(overall, {x: \"dist\", y: \"meanScore\", fill: d => d.dist == \"kernel\" ? \"blue\" : \"grey\"}),\n Plot.ruleY([overall[0].meanScore], {stroke: \"blue\", strokeDasharray: \"2,2\"}),\n Plot.ruleY([0])\n ]\n })\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":87,"value":"### Comparing Forecast Horizons","pinned":false,"hidden":false,"mode":"md"},{"id":92,"value":"horizonResults = {\n const grouped = d3.flatGroup(scoreResults, d => d.tp)\n return grouped.map(d => {\n const meanDistScore = d3.mean(d[1], d => d.distScore)\n const meanNormalScore = d3.mean(d[1], d => d.normalScore)\n return {tp: d[0], meanDistScore, meanNormalScore, diff: meanDistScore-meanNormalScore}\n })\n return grouped\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":526,"value":"{\n const tpPlot = Plot.plot({\n height: 260,\n width: width/2,\n marks: [\n Plot.lineY(horizonResults, {x: \"tp\", y: \"meanNormalScore\", stroke: \"grey\", strokeDasharray: \"1,3\"}),\n Plot.lineY(horizonResults, {x: \"tp\", y: \"meanDistScore\", stroke: \"blue\", strokeDasharray: \"1,3\"}),\n Plot.dot(horizonResults, {x: \"tp\", y: \"meanNormalScore\", stroke: \"grey\"}),\n Plot.dot(horizonResults, {x: \"tp\", y: \"meanDistScore\", stroke: \"blue\"}),\n ]\n })\n\n const tpPlot2 = Inputs.table(horizonResults, {\n width: width/2-100, height: 260, format: {\n diff: d => htl.html`

${d.toFixed(3)}

`\n }})\n \n return plotFlex([tpPlot, tpPlot2])\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":312,"value":"### Over Time","pinned":false,"hidden":false,"mode":"md"},{"id":251,"value":"viewof plotTp = Inputs.range([0, tp], {label: \"Plot tp\", step: 1})","pinned":false,"hidden":false,"mode":"ojs"},{"id":321,"value":"{\n const resultsTp = scoreResults.filter(d => d.tp == plotTp)\n return Plot.plot({\n height: 300,\n marginLeft: 70,\n marks: [\n Plot.areaY(resultsTp, {x: \"baseT\", y: \"normalScore\", fill: \"rgba(128,128,128,.2)\"}),\n Plot.areaY(resultsTp, {x: \"baseT\", y: \"distScore\", fill: \"rgba(0,0,255,.2)\"}),\n Plot.ruleY([0]),\n Plot.ruleX([plotT])\n ]\n })\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":249,"value":"### Specific forecast","pinned":false,"hidden":false,"mode":"md"},{"id":253,"value":"viewof plotT = {\n let extent = d3.extent(scoreResults, d => d.baseT)\n extent[1] = extent[1] - plotTp \n return Inputs.range(extent, {label: \"Plot t\", step: 1})\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":263,"value":"plotForecast = distForecasts.find(d => d.baseT == plotT && d.tp == plotTp)","pinned":false,"hidden":false,"mode":"ojs"},{"id":258,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":395,"value":"{\n const res = scoreResults.find(d => d.baseT == plotT && d.tp == plotTp)\n const tableResults = [{normalScore: res.normalScore, distScore: res.distScore, diff: res.distScore - res.normalScore}]\n return Inputs.table(tableResults, {\n width: 360,\n format: {\n diff: d => htl.html`

${d.toFixed(3)}

`\n }\n })\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":594,"value":"## Point Peak Scores\n\nA key idea behind EpiForecast is that the user can use their domain knowledge to choose appropriate information from the interpretable forecasts to get a better intuition of the potential futures. One intuitive reading of the shaded forecast area is that each peak represents the center of a potential future identified by the method. Here we compare the best case scenario (the user always picks the best potential future), the worst case scenario (the user always picks the worst potential future), and the simplex point forecast. To do this we treat each peak as a point forecast and compare the mean absolute error (MAE) between the best performing peak, the worst performing peak, and the point forecast generated using the simplex method. ","pinned":false,"hidden":false,"mode":"md"},{"id":598,"value":"horizonPeakResults = {\n const grouped = d3.flatGroup(scoreResults, d => d.tp)\n return grouped.map(d => {\n const meanPointError = d3.mean(d[1], d => d.pointError)\n const meanDistPeakErrorBest = d3.mean(d[1], d => d.distPeakErrorBest)\n const meanDistPeakErrorWorst = d3.mean(d[1], d => d.distPeakErrorWorst)\n const meanDistPeakErrorMean = d3.mean(d[1], d => d.distPeakErrorMean)\n //const meanMeanPeakError\n return {tp: d[0], meanPointError, meanDistPeakErrorBest, meanDistPeakErrorWorst, meanDistPeakErrorMean}\n })\n return grouped\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":607,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":59,"value":"## Libraries and Helper Functions","pinned":false,"hidden":false,"mode":"md"},{"id":572,"value":"function getPeaks(pdf, domain, n=100) {\n const points = funcPoints(pdf, domain, n)\n\n const peaks = []\n for (let i = 1; i < points.length-1; i++) {\n if (points[i][1] > points[i-1][1] && points[i][1] > points[i+1][1]) {\n peaks.push(points[i])\n }\n }\n return peaks\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":550,"value":"function plotFlex(plots) {\n const div = htl.html`

`\n div.style = \"display: flex; justify-content: space-around\"\n plots.forEach(plot => div.appendChild(plot))\n return div\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":433,"value":"function naiveNegativeLinearScore(dist, x) {\n return -dist(x)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"value":"function negativeLogScore(dist, x, minPd=1e-07) {\n return -Math.log(Math.max(dist(x), minPd))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":154,"value":"function funcPoints(f, extent, n=100) { \n const step = (extent[1]-extent[0])/n\n const points = []\n for (let i = 0; i < n; i++) {\n const val = extent[0] + i*step\n points.push([val, f(val)])\n }\n return points\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":107,"value":"// TODO: Add to forecast library, then remove\nfunction kdeFunctions(XW, opts={}) {\n opts = Object.assign({kernelWidth: null, n: 80, c: 0.0000001}, opts)\n let {kernelWidth, n, c} = opts\n\n if (kernelWidth == null) {\n kernelWidth = Math.max(d3.deviation(XW, d => d[0]), 0.000001)*2\n }\n\n const norm = gaussian(0, 1)\n\n const s = 1 / (norm.pdf(0)*XW.length)\n const h = Math.sqrt((-2/kernelWidth**2)*Math.log(c*Math.sqrt(2*Math.PI)))\n\n const kernelPdfFunctions = XW.map(([x,w]) => {\n return d => s * w * norm.pdf(h*(d-x))\n })\n\n const vRange = d3.extent(XW, d => d[0])\n const domain = [vRange[0]-kernelWidth, vRange[1]+kernelWidth]\n\n const pdfFunction = d => d3.sum(kernelPdfFunctions, f => f(d))\n return {pdf: pdfFunction, kernelPdfs: kernelPdfFunctions, relDomain: domain}\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":96,"value":"function pdf(dist, v) {\n const i = dist.findIndex(d => d.v > v)\n if (i == 0 || i == -1) return 0\n \n const window = [dist[i-1], dist[i]]\n const p = (v - window[0].v)/(window[1].v-window[0].v)\n return window[0].p + p*(window[1].p-window[0].p)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":177,"value":"function areaUnderFunc(f, domain) {\n const data = funcPoints(f, domain)\n let totalArea = 0\n const xf = d => d[0]\n const yf = d => d[1]\n for (let i = 0; i < data.length-1; i++) {\n const h = xf(data[i+1])-xf(data[i])\n const a = yf(data[i])\n const b = yf(data[i+1])\n totalArea += h*(a + b)/2\n }\n return totalArea\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":78,"value":"function areaUnder(data, xf, yf) {\n let totalArea = 0\n for (let i = 0; i < data.length-1; i++) {\n const h = xf(data[i+1])-xf(data[i])\n const a = yf(data[i])\n const b = yf(data[i+1])\n totalArea += h*(a + b)/2\n }\n return totalArea\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":570,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":81,"value":"gaussian = (await import('https://cdn.skypack.dev/gaussian@1.3.0?min')).default","pinned":false,"hidden":false,"mode":"ojs"},{"id":62,"value":"fc = import(`https://episphere.github.io/forecast/src/index.js`)","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}