Don't forget snap, crackle, and pop! I design machinery for entertainment, and we deal with jerk, the third derivative of displacement, and sometime snap, the fourth derivative. I never thought of them applied this way, but I guess it makes sense that you can now that I think about it. Interesting.
I just started using R/RStudio in the past few months. I'd be very interested in the J and R connection. Thank you! On Sat, 28 Mar 2020 at 07:18, Dimitri Georganas <[email protected]> wrote: > 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 > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
