Hi, I've been modeling the spread of the disease using a Gompertz function and its 1st, 2nd and 3d derivative to determine rate, acceleration and 'jerk', but I couldn't get the R addon to work with my J installation, so I did this part in R. (I couldn't find Levenberg Marquardt curve fitting libraries/addons in J) Then I used the a,b,c I got from the R code in a J script to estimate the number of IC beds required, assuming a 10% critical cases scenario. Well, Houston, we have a problem.
I could port the R code to J if I could either get the R addon to work or find a Levenberg-Marquardt implementation in J. IC impact code https://github.com/dimgeo/IC-covid R code with Gompertz curve fitting and derivative https://github.com/dimgeo/covid Best regards, Dimitri On Fri, Mar 27, 2020 at 7:24 PM Robert Herman <[email protected]> wrote: > I found this article: > > Modeling the COVID-19 Outbreak with J > < > https://datakinds.github.io/2020/03/15/modeling-the-coronavirus-outbreak-with-j > > > > I took the author's code, which is a modified SEIR model (The_SEIR_model > < > https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model > >), > and updated some of the data and recalculated the variables. I then > pilfered some of Michal's (Tangentstorm, you rock!) code for animating > things in J from his talks. This is helping me to push my learning of J, > and keep it interesting. > > Two disclaimers: The author of the linked article doesn't claim to be an > epidemiologist, and I am only doing this to learn J, and cannot comment on > any of the validity of the algorithm or simulations it produces. > > Now, on to my two issues (full code at end): > > 1. I know I have to figure out a way to have the* infect_display^:144 > starting population* call feed the *update* loop, otherwise it is trying to > run all 144 days of simulation in one call to update at a frequency set by > *wd > 'timer 600'*. I am working on this, and realize I somehow have to run it > outside of update, and then feed update "a day at a time". The goal is to > see the same person (grid location) evolve over the simulation. I know this > is simple, but i haven't gotten there yet. Any tips will be appreciated! > > 2. When I had initially set up the code to run 25 people instead of 144, > the color index I had setup with the DSEIR variable exhibit two strange > glitches. First, squares that were supposed to be 1=White, were rendering > as 0=Black=Dead, and sometimes colors appeared "blended". An yellowish, > off-white for yellow, and at times bright yellow. When I went to 144 people > and it is trying to run in update, the grid appears to be correctly > indexing and not moving askew from frame to frame, but the colors are > inconsistent. > > BTW, thanks to Linda's article in Journal of J - Vol.4, No.1 for the RGB > list to then pull my desired colors from in the DSEIR variable (Dead, > Susceptible, Exposed, Infectious, Recovered): > > *RGB=: (#:i.8) { 0 255* > > > *DSEIR =: (0 7 6 4 2 { RGB)* > > Here is my integration of the two sources of code to try and animate the > outbreak per the article, and thanks!: > > *load 'viewmat'* > > *coinsert'jgl2'* > > > *wd 'pc Covid-19-Sim closeok' NB. parent control (window) named 'w0'* > > *wd 'minwh 500 500; cc g0 isidraw;' NB. add an 'isidraw' child control > named 'g0'* > > *wd 'pshow; pmove 40 510 0 0' NB. show the window at the given > coordinates.* > > *wd 'sm focus term' NB. session manager: bring terminal to front* > > *wd 'psel Covid-19-Sim; ptop' NB. bring our window to front* > > > *vmcc =: viewmatcc_jviewmat_ NB. viewmat to a child control* > > > *step =: render @ update NB. each step, we'll call those two in sequence* > > > *RGB=: (#:i.8) { 0 255* > > > *DSEIR =: (0 7 6 4 2 { RGB)* > > > *update =: verb define* > > * im =: 12 12 $ infect_display^:144 starting_pop* > > *)* > > > *render =: verb define* > > * DSEIR vmcc im;'g0'* > > * glpaint''* > > *)* > > > *sys_timer_z_ =: step_base_* > > *wd 'timer 600'* > > > *ppl =: 144 NB. How many people are in our simulation?* > > > *Beta =: 0.106271 NB. The rate at which susceptible people become exposed. > RPH - dependent on R0 = 2.2 3/25/2020* > > *Sigma =: 0.066967 NB. The rate at which exposed people become infectious. > RPH - (1-sigma)^10=0.5, 10 is incubation period.* > > *Gamma =: 0.056126 NB. The rate at which infectious people recover. RPH - > adjusted to 12 days, so (1-Gamma)^12=0.5.* > > *Xi =: 0.02284 NB. The rate at which recovered people lose immunity and > become susceptible. RPH - 30 days* > > * NB. is partially based upon man on Diamond Princess who presented > symptoms around 30 days otherwise made up.* > > *R =: 2.2 NB. How many people will 1 person infect? RPH - based on latest > CDC estimate 3/25/2020* > > *CFR =: 0.015 NB. Case Fatality Rate - RPH based upon 0.25% to 3% from CDC, > so 1.5% is higher than average.* > > *D =: 0.0011619 NB. Death rate per day while symptomatic - RPH based upon > 13 days to death as seen and per CDC.* > > > *ClosenessStdDev =: 5 % 3* > > *ClosenessMean =: 1* > > > *rand_norm =: 3 : '(2&o. 2p1 * ? y) * %: _2 * ^. ? y'* > > *closeness =: 4 %~ 2 %~ (+ |:) ClosenessMean + ClosenessStdDev * rand_norm > (ppl,ppl) $ 0* > > *risk =: closeness - closeness * =i.ppl * > > *NB. risk =: (ppl,ppl)$0* > > > *NB. 0 = Dead* > > *NB. 1 = Susceptible* > > *NB. 2 = Exposed* > > *NB. 3 = Infectious* > > *NB. 4 = Recovered* > > *starting_pop =: 2 , (ppl - 1) $ 1* > > > *infect =: 3 : 0* > > * can_spread =. (1&< *. 4&>) y* > > * susceptible =. y = 1* > > * exposed =. y = 2* > > * infectious =. y = 3* > > * recovered =. y = 4* > > > * infectiousness =. can_spread ,./ . * risk* > > > * susceptible_to_exposed =. (Beta * susceptible) * (+/ can_spread) %~ +/ |: > susceptible * infectiousness* > > * exposed_to_infectious =. Sigma * exposed* > > * infectious_to_recovered =. Gamma * infectious* > > * infectious_to_dead =. D * infectious* > > * recovered_to_susceptible =. Xi * recovered* > > * 1 (I.recovered_to_susceptible>?ppl$0)} 4 > (I.infectious_to_recovered>?ppl$0)} 0 (I.infectious_to_dead>?ppl$0)} 3 > (I.exposed_to_infectious>?ppl$0)} 2 (I.susceptible_to_exposed>?ppl$0)} y* > > *)* > > > *infect_display =: 3 : 0* > > * out =. infect y* > > * smoutput out* > > * smoutput > > ('Susceptible';'Exposed';'Infectious';'Recovered';'Dead'),:(+/1=out);(+/2=out);(+/3=out);(+/4=out);(+/0=out)* > > * out* > > *)* > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
