31:T458,md`# "Live" Logistic Coronavirus Death Counter This page gives some representation of how often deaths are occuring due to Coronavirus. Most [live counters](https://www.worldometers.info/) of events just show a counter that increments at a regular rate, e.g. once every 15 seconds. Naturally occuring rare events instead show [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution) like statistics, with [exponential](https://en.wikipedia.org/wiki/Exponential_distribution) waiting times. This page shows a "live" counter that estimates what the current number of deaths would be from an logistically growing process. For more about logistic growth, see [this notebook](https://observablehq.com/@alemi/logistic-growth-plots). It uses the latest data from [Johns Hopkins](https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6), and a simple model. Below you can change the country or seed to vary the results slightly. I emphasize that this isn't actually a live number, but gives a sense of how the deaths would be spread or occuring throughout the day. `32:T45a,md`The current average time between deaths is: ${formatTimeDelta( avgTime )}. This means there are on average ${d3.format(".2f")( (1000. * 60 * 60) / avgTime )} deaths each hour. ${ country == "US" ? md`To put that in perspective in [2017](https://www.cdc.gov/nchs/fastats/deaths.htm) on average ${d3.format( ".2f" )( 2813503 / (365 * 24) )} Americans died per hour from any cause. The leading cause of death was heart disease from which ${d3.format( ".2f" )(647457 / (365 * 24))} people died each hour.` : `` } The latest estimate for the base rate is ${d3.format(".3f")( currentAlpha )} / person / day meaning the doubling time is currently ${d3.format(".3f")( Math.log(2) / (currentAlpha * (1 - currentDeaths / currentCap)) )} days, and the estimate for the carrying capacity is ${d3.format(".0f")( currentCap )}. The estimated peak rate is therefore ${d3.format(".2f")( (currentAlpha * currentCap) / 4 / 24 )} per hour, which should ${ estimatedPeakDate >= currentDate ? "occur" : "have occured" } at ${estimatedPeakDate.toLocaleString()} `33:T948,md`# Model I wanted to have a "live" ticker that more accurately represented the spacing between events than just a regularly ticking clock. As our base model, we will assume the Coronavirus deaths grow logistically. This means we will assume that the probability of a death in a short amount of time ${tex`\delta t`} is proportional to some base rate per person ${tex`\lambda`} times the number of current deaths (${tex`N`}) times a factor of how close we are the the *carrying capacity* (${tex`(1-N/K)`}): ${tex.block`P(\text{death}) \sim \lambda N \left(1 - \frac{N}{K}\right) \delta t`} What do we set ${tex`\lambda`} and ${tex`K`} to? We want our model to be accurate, so it would be nice if existing data could inform our decision. Additionally, I want this to be like a true counter in the sense that it only ever can increase, never decrease. We will use the data from [Johns Hopkins Github](https://github.com/CSSEGISandData/COVID-19) which updates at midnight UTC every night. Broadly we are going to want to: every night, look at the data, form some estimate of how many deaths we expect over the next 24 hours, and in concert with our current value, adjust our rate to approximately hit that target. To do this we need to characterize how many deaths we expect to happen in any given amount of time. What we want is the expected number of deaths at any given time. Solving the above differential equation we get for the expected behavior: ${tex.block`N(t) = \frac{K}{1 + e^{-\lambda(t - t_0)}}`} How do we estimate the parameters? We'll notice that if we rearrange the original differential equation we can obtain: ${tex.block`\frac{\Delta N}{N} = \frac{\lambda}{K} ( K - N )`} This form is useful because it tells that we expect a *linear* relationship between the fractional change in the number of deaths in a period, and the number of deaths. This allows us to read off the paramters. For this notebook, we'll take these estimate from successive days, where the deltas are measured over the course of a week. If we want to express our rate in units of days, this means we have the following system of equations to determind our parameters. ${tex.block`\frac{N(t) - N(t-7)}{7 N} = \frac{\lambda}{K} (K - N(t)) \\ \frac{N(t-1) - N(t-8)}{7 N} = \frac{\lambda}{K}(K - N(t-1)) `} which can be solved to estimate the parameters at any moment. `34:T52d,function statsRecent(time) { let recents = deathCount.filter(x => x.time <= time).reverse(); if (recents.length == 0) { return { date: undefined, time: undefined, dateStr: undefined, deaths: 0, xn: 0, xnm1: undefined, xnm7: undefined, xnm8: undefined, alpha: priorRate, k: priorCount }; } if (recents.length < 9) { return { date: recents[0].date, time: recents[0].time, dateStr: recents[0].dateStr, deaths: recents[0].deaths, xn: recents[0].deaths, xnm1: undefined, xnm7: undefined, xnm8: undefined, alpha: priorRate, k: priorCount }; } let r = { date: recents[0].date, time: recents[0].time, dateStr: recents[0].dateStr, deaths: recents[0].deaths, xn: recents[0].deaths, xnm1: recents[1].deaths, xnm7: recents[7].deaths, xnm8: recents[8].deaths }; r.alpha = (1. / 7) * (1 + (-r.xnm7 * r.xnm1 * r.xnm1 + r.xnm8 * r.xn * r.xn) / (r.xnm1 * (r.xnm1 - r.xn) * r.xn)); r.k = (r.xnm7 * r.xnm1 * r.xnm1 - r.xnm1 * r.xnm1 * r.xn - r.xnm8 * r.xn * r.xn + r.xnm1 * r.xn * r.xn) / (r.xnm7 * r.xnm1 - r.xnm8 * r.xn); // r.k = (r.alpha * r.xn * r.xn) / (r.xnm7 - r.xn + r.alpha * r.xn); return r; }35:T5bb,events = { let arr = new Array(); let time = deathCount[0].time; const startTime = time; let n = deathCount[0].deaths; const random = d3.randomExponential.source(seedrandom(seed))(1.0); const lastTime = deathCount.last().time; const finalTime = deathCount.last().time + 14 * millisecondsPerDay; let days = 0; arr.push({ time: time, deaths: n, date: new Date(time), deltaDays: (time - startTime) / millisecondsPerDay }); while (time < finalTime) { const stats = statsRecent(time); const alpha = stats.alpha; const k = stats.k; let rate = Math.max(0, alpha * n * (1 - n / k)); const dayEnd = startTime + (days + 1) * millisecondsPerDay; let dt = (dayEnd - time) / millisecondsPerDay; dt = 1; let n1 = stats.deaths; let n2 = n; let fudgeAlpha = time < lastTime ? Math.log((n1 * (k - n2)) / ((k - n1) * n2)) / dt : 0; while (time < dayEnd) { let currentRate = Math.max(0, (alpha + fudgeAlpha) * n * (1 - n / k)); let waitTime = (random() / currentRate) * millisecondsPerDay; if (time + waitTime <= dayEnd) { time = time + waitTime; n = n + 1; arr.push({ time: time, deaths: n, date: new Date(time), deltaDays: (time - startTime) / millisecondsPerDay }); //yield arr; } else { time = dayEnd; days = days + 1; break; } } } return arr; }36:[{"id":"6e4514e5067a7bfd7c568b6e6978f5c099aee742a27c59604c5454fba67bf573ab19dc0c4e05465d2c23265fe8281fffd9dd670a11ccdfcf9c156b0295686508","href":"https://static.observableusercontent.com/files/6e4514e5067a7bfd7c568b6e6978f5c099aee742a27c59604c5454fba67bf573ab19dc0c4e05465d2c23265fe8281fffd9dd670a11ccdfcf9c156b0295686508","name":"378799__kgeshev__bell.wav","type":"audio/vnd.wave","lastModified":1585344298072.403,"size":1320966}] 2f:{"title":"\"Live\" Logistic Coronavirus Death Counter","cells":[{"id":880,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":713,"value":"md`\nThe current counter of ${country} Corona deaths is ${nextEvent.deaths - 1},\nwith the next death expected @ ${nextEvent.date} in ${formatTimeDelta(\n remainingTime\n)}.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":729,"value":"html`

${currentDeaths} Deaths

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":819,"value":"html`

Next death in ${formatTimeDelta(\n remainingTime\n)}

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":988,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1159,"value":"viewof playSound = checkbox({\n description: \"Enable Bell Tones\",\n options: [{ value: \"toggle\", label: \"Play Sound\" }],\n value: \"toggle\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":40,"value":"viewof country = autoSelect({\n options: countries,\n value: \"US\",\n title: \"Country\",\n description: \"Search for a country.\",\n placeholder: \"Search for a country . . .\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":822,"value":"html`

Recent Events

\n\n\n

Death	Date\n\n${events\n .slice(nextIndex - 15, nextIndex)\n .reverse()\n .map(\n x =>\n \"
\" + x.deaths + \"	\" + d3.timeFormat(\"%m-%d @ %H:%M:%S\")(x.date)\n )}\n

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":927,"value":"md`Below is a plot of how accurate this approach has been over time. Blue is actual data, red is the state of the ticker.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":174,"value":"{\n const obs = vl\n .markLine({ color: \"steelblue\", strokeWidth: 7 })\n .data(deathCount)\n .title(\"Deaths\")\n .encode(\n vl.x().fieldT(\"date\"),\n vl.y().fieldQ(\"deaths\")\n // .scale({ type: \"log\" })\n );\n\n const guess = vl\n .markLine({ color: \"firebrick\", interpolate: \"step-after\" })\n .data(events)\n .title(\"Deaths\")\n .encode(\n vl.x().fieldT(\"date\"),\n vl\n .y()\n .fieldQ(\"deaths\")\n .scale({ type: \"log\" })\n );\n\n return vl.layer(obs, guess).render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1264,"value":"md`Changing the seed changes the curve slightly, any string will do.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":282,"value":"seed = \"53090\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1273,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":1291,"value":"md`\nIn order to keep the model in sync with the actual observations, we also adjust the expected rate to ensure that both the current state of our model as well as the latest measurements of the data should agree at the end of the day.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1207,"value":"md`## Code\n\nImplementation, not terribly organized.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"deaths = d3.csv(\n \"https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv\"\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"dates = {\n let s = new Set();\n for (var k of deaths) {\n for (var v in k) {\n s.add(v);\n }\n }\n for (var v of [\"Province/State\", \"Country/Region\", \"Lat\", \"Long\"]) {\n s.delete(v);\n }\n return Array.from(s);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1220,"value":"countries = {\n let s = new Set();\n for (var k of deaths) {\n s.add(k[\"Country/Region\"]);\n }\n return Array.from(s);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":233,"value":"deathCount = dates\n .map(date => ({\n dateStr: date,\n date: dateStrToDate(date).addDays(1),\n time: dateStrToDate(date)\n .addDays(1)\n .getTime(),\n deaths: deaths\n .filter(x => x[\"Country/Region\"] == country)\n .map(x => parseInt(x[date]))\n .reduce((acc, v) => v + acc)\n }))\n .drop(x => x.deaths > 0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":765,"value":"currentDate = new Date(currentTime)","pinned":false,"hidden":false,"mode":"ojs"},{"id":255,"value":"currentTime = {\n var eventPK = bisector(events, Date.now());\n while (true) {\n let now = Date.now();\n if (playSound && now > events[eventPK].time) {\n audio.play();\n ++eventPK;\n }\n yield Promises.tick(1000 / 60, Date.now());\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":285,"value":"function statsLess(time) {\n return deathCount.filter(x => x.time <= time).last();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1466,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":246,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":689,"value":"bisector = d3.bisector(function(d) {\n return d.date;\n}).right","pinned":false,"hidden":false,"mode":"ojs"},{"id":695,"value":"nextEvent = events[nextIndex]","pinned":false,"hidden":false,"mode":"ojs"},{"id":748,"value":"nextIndex = bisector(events, currentTime)","pinned":false,"hidden":false,"mode":"ojs"},{"id":991,"value":"currentDeaths = nextEvent.deaths - 1","pinned":false,"hidden":false,"mode":"ojs"},{"id":490,"value":"remainingTime = nextEvent.time - currentTime","pinned":false,"hidden":false,"mode":"ojs"},{"id":1570,"value":"estimatedPeak = (Math.log(currentCap - 1) +\n Math.log(\n (currentCap - deathCount.last().deaths) /\n ((currentCap - 1) * deathCount.last().deaths)\n )) /\n currentAlpha","pinned":false,"hidden":false,"mode":"ojs"},{"id":1580,"value":"estimatedPeakDate = new Date(\n deathCount.last().time + estimatedPeak * millisecondsPerDay\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1567,"value":"estimateEnd = (Math.log(currentCap - 1) +\n Math.log(\n (99 * (currentCap - deathCount.last().deaths)) /\n ((currentCap - 1) * deathCount.last().deaths)\n )) /\n currentAlpha","pinned":false,"hidden":false,"mode":"ojs"},{"id":187,"value":"priorCount = 1000000","pinned":true,"hidden":false,"mode":"ojs"},{"id":189,"value":"priorRate = 0.9 / 7","pinned":true,"hidden":false,"mode":"ojs"},{"id":1462,"value":"currentAlpha = statsRecent(currentTime).alpha","pinned":false,"hidden":false,"mode":"ojs"},{"id":1461,"value":"currentCap = statsRecent(currentTime).k","pinned":false,"hidden":false,"mode":"ojs"},{"id":977,"value":"currentActiveRate = currentAlpha *\n currentDeaths *\n (1 - currentDeaths / currentCap)","pinned":false,"hidden":false,"mode":"ojs"},{"id":983,"value":"avgTime = 1. / (currentActiveRate / millisecondsPerDay)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1564,"value":"{\n const obs = vl\n .markLine({ color: \"steelblue\", strokeWidth: 5 })\n .data(deathCount)\n .title(\"Deaths\")\n .encode(vl.x().fieldT(\"date\"), vl.y().fieldQ(\"deaths\"));\n\n const guess = vl\n .markLine({ color: \"firebrick\", interpolate: \"step-after\" })\n .data(events)\n .title(\"Deaths\")\n .encode(vl.x().fieldT(\"date\"), vl.y().fieldQ(\"deaths\"));\n\n return vl.layer(obs, guess).render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":336,"value":"md`## Library\n\nUtility functions.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":972,"value":"millisecondsPerDay = 24 * 60 * 60 * 1000","pinned":false,"hidden":false,"mode":"ojs"},{"id":619,"value":"(Date.prototype.addDays = function(days) {\n var date = new Date(this.valueOf());\n date.setDate(date.getDate() + days);\n return date;\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":339,"value":"(Array.prototype.last = function() {\n return this[this.length - 1];\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":347,"value":"(Array.prototype.drop = function(cond) {\n return this.slice(this.findIndex(cond));\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":411,"value":"dateStrToDate = d3.utcParse(\"%m/%d/%y\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":795,"value":"formatTimeDelta = function(diff) {\n var diff = diff;\n var output = \"\";\n if (diff > 24 * 60 * 60 * 1000) {\n let days = Math.floor(diff / (24 * 60 * 60 * 1000));\n output += days + \" days\";\n diff = diff - days * 24 * 60 * 60 * 1000;\n }\n if (diff > 60 * 60 * 1000) {\n let hours = Math.floor(diff / (60 * 60 * 1000));\n output += (output.length > 0 ? \" \" : \"\") + hours + \" hours\";\n diff = diff - hours * 60 * 60 * 1000;\n }\n if (diff > 60 * 1000) {\n let minutes = Math.floor(diff / (60 * 1000));\n output += (output.length > 0 ? \" \" : \"\") + minutes + \" minutes\";\n diff = diff - minutes * 60 * 1000;\n }\n output +=\n (output.length > 0 ? \" \" : \"\") + d3.format(\".2f\")(diff / 1000) + \" seconds\";\n return output;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1078,"value":"audio = html`