El factor de crecimiento siempre tiende a uno en una función polinómica

A ver, si tienes una función que crece como \(f(x) = b x^n\), y te pones a calcular el factor de incremento \(F(x) = f(x+a)/f(x)\), es obvio que tiende a uno:

$$F(x \to \infty) = (1+a/x)^n$$

Por supuesto, cualquier regimen exponencial \(f(x) = b e^{k x^n}\) tiene un factor de incremento \(F(x) = e^{k [(x+a)^n – x^n]}\) que se va a infinito para \(n>1\) mientras que para \(n=1\) tiende a $$F(x \to \infty) =e^{k a}$$ un límite que, al ser el exponente positivo, es siempre mayor que uno.

Así, la convergencia a un factor de crecimiento unidad tan solo nos indica que hemos abandonado el régimen exponencial. La “fase de ralentización”. De hecho, ni siquiera tenemos garantizado haber entrado en el regimen polinomico, podriamos estar en el subexponencial, n <1 del caso anterior. Por ejemplo exp(sqrt(x)).

Extensiones del modelito no compartimental de ayer

Primero un par de avisos de interes general:

Mirando esos apuntes y modelos se ve que los modelos compartimentales tienen muchas posibilidades, especialmente cuando comienzas a permitir que una valvula tenga dependencia funcional arbitraria, y no simplemente productos.

Pero bueno, ya que estaba con la idea no compartimental de ayer, voy a dejar aun anotadas dos o tres cosas que se me quedaron en el cuaderno.

La idea de poner un tiempo de incubacion gaussiano o con concentracion de delta de dirac no parece contraria a la intuicion, pero hacer eso con el tiempo de resolucion/removal no parece tan logico, habria que ir poniendo otras distribuciones d(t), Weibulls, cosas asi. Ahora, a falta de otro criterio conviene que podamos hacer analiticamente la integral, asi que vamos a hacer algo que suena cutre pero no deja de ser parecido a lo de los compartimentos: asumir que la probabilidad de detectar al infectado funciona con una distribucion exponencial. Me direis que para ese viaje no hacian falta alforjas, pero es bastante curioso porque eso nos lleva a

$$
(1 – \int_0^\infty d(t) e^{-rt} dt) = 1 – \int_0^\infty d(t) e^{-rt} dt = { r \over \lambda +r}
$$

con \(\lambda = 1/\tau_d\), y este seria el tiempo medio de deteccion, asi que $$R0 = k \tau_d$$.

La ecuacion de ayer para relacionar k con el exponente de crecimiento inicial se simplifica casi por sorpresa:
$$ k { r \over \lambda +r} = r e^{r t_l} $$
y ahora la solucion grafica es encajar dos una hiperbola \({ k \over \lambda +r}\) con la exponencial pura \(e^{r t_l}\) que condensa el tiempo de latencia. Este segundo termino vale siempre 1 en el origen, asi que la condicion de existencia de solucion es que la hiperbola debe ser mayor que uno en el origen, usease
$$ {k \over \lambda} > 1$$
¡que vuelve a ser la condición basica de la teoria epidemiologica: R0 > 1! Al menos esto parece que nos da bien.

Si no tuvieramos tiempo de incubacion la solucion de hecho seria bastante trivial:

$$r = k – {1 \over \tau_d}$$

Indicando que el exponente que vemos en las gráficas que salen por internet hay que esperar que tenga al menos dos contribuciones: la de la velocidad de contagio y la del tiempo de detección. Y por eso es por lo que tenemos siempre el debate sobre si queremos meterle presión a bajar una cosa o la otra.

Otro apunte. Si lo que queremos es precisamente ver la sensibilidad de r a cada uno de estos parametros puede que haya casos en los que en vez de resolver explicitamente (usando más aproximaciones) nos podamos apañar con los teoremas de derivadas parciales de funciones implicitas, usease consideramos

$$
0 = F(r(k,\lambda),k,\lambda) = {k \over \lambda +r} e^{rt_l}
$$

y aplicamos las reglas de encadenado aquellas de segundo de carrera

$$
{\partial r \over \partial k} = – {{\partial F \over \partial k} \over {\partial F \over \partial r} }
$$

etcetera… si conseguimos recordarlas y asegurarnos de que se cumplen las condiciones 😀

Mi hot take de la expansion exponencial del coronavirus

Voy a poner por aqui mi variante de SIR con tiempos de incubacion y tiempos de deteccion (respuesta, resolucion) hasta que encuentre algun modelo oficial que citar.

Para tener expansion exponencial, la masa susceptible es infinita, asi que en realidad lo que nos interesa es modelar Expuestos — > Infecciosos — > Resueltos. Si fuera un modelo compartimental seria un caso extremo de SEIR, como me apuntaba Francis en twitter.

Si en vez de usar grifos, o como se llamen las valvulas entre compartimentos, modelamos el salto E–>I con una funcion de latencia \(l(t)\) y el salto I–>R con una funcion de deteccion \(d(t)\), creo que el incremento de infectados seria algo asi como:

$$
\delta I_+(t)=\int_0^\infty l(t’) k I(t-t’) dt’
$$

Donde el contagio propiamente va con una velocidad k tal que
$$
\delta E = k I
$$

Con lo que tengo mas dudas es con la deteccion, porque hay que recordar de alguna manera en que momento entro cada paciente en la caja de infectados. A vuelapluma se me ocurre que seria algo asi:

$$
\delta I_-(t) = \int_0^\infty \int_0^\infty d(t”) l(t’) k I(t-(t’+t”)) dt’ dt”
$$

El total de casos, C(t), seria la suma de tan solo la parte positiva de la aportacion, mientras que la equacion diferencial para los infectados seria la union de los dos aportes

$$
I'(t) = \delta I_+(t) – \delta I_-(t)
$$

Si asumimos que la solucion es exponencial, \(I(t) = I_0 e^{r t}\), entonces la “convolucion” del termino de resolucion es separable:
$$
\delta I_-(t) = \int_0^\infty \int_0^\infty d(t”) l(t’) k I_0 e^{rt} e^{-rt’} e^{-rt”} dt’ dt”
$$
y se simplifica todo bastante
$$
I_0 r e^{r t} = \int_0^\infty l(t’) k I_0 e^{rt} e^{-rt’} dt’- \int_0^\infty \int_0^\infty d(t”) l(t’) k I_0 e^{rt} e^{-rt’} e^{-rt”} dt’ dt”
$$

$$
r = \int_0^\infty l(t’) k e^{-rt’} dt’- \int_0^\infty \int_0^\infty d(t”) l(t’) k e^{-rt’} e^{-rt”} dt’ dt”
$$

$$
r = k \big( \int_0^\infty l(t’) e^{-rt’} dt’ \big) \big(1 – \int_0^\infty d(t”) e^{-rt”} dt”\big)
$$

Por otro lado el dato medido en prensa son los casos totales, vendria a ser

$$
C(t) = {k\over r} I_0 \big( \int_0^\infty l(t’) e^{-rt’} dt’ \big) e^{rt}
$$

Si vamos al caso particular mas sencillo en el que las funciones l(t) y d(t) se concentran simplemente en unos puntos t_l y t_d, las integrales se resuelven con la delta de Dirac y nos queda

$$
r = k e^{-rt_l} \big(1 – e^{-rt_d}\big)
$$

$$
C(t) = {k\over r} I_0 e^{-rt_l} e^{rt} = {k\over r} I_0 e^{r (t-t_l)}
$$

lo que parece razonable teoricamente, porque al menos se cumple la condicion universal de que para que haya solucion tiene que ocurrir que el termino R0 sea mayor que uno. Explicitamente, si queremos resolver para encontrar \(r\) ponemos la ecuacion separada en dos terminos
$$
k (1 – e^{-r t_d}) = r e^{r t_l}
$$
tenemos dos funciones sin puntos de inflexion y distinta concavidad, asi que para que se crucen necesitamos que la derivada en en el origen del termino izquierdo sea mayor que la del termino derecho. La de este ultimo es siempre uno, independientemente del tiempo de latencia de la enfermedad, mientras que la del termino izquierdo es \(k t_d\), asi que existe solucion exponencial solo cuando se cumple que
$$
k t_d > 1
$$
Y esto no es otro que el parametro R0: la velocidad de contagio multiplicado por el tiempo que se tarda hasta que se detecta y resuelve.

En el caso real conocemos \(r \approx 0.33\) pero no tenemos ni idea de los otros tres parametros, k, t_l y t_d.

NOTA:
Si asumimos que R0 siempre crece con el tiempo de deteccion, da la impresion de que en este modelo, y posiblemente de modo mas general, dado un valor para el tiempo de latencia podemos dar una cota inferior para R0, que sería \(e^{r t_l}\). Pero para r=0.34 el resultado no es muy consistente, una latencia de cinco dias corresponderia ya a una cota R0 > 5.47; para obtener R0 razonable entre uno y 2.5 habria que asumir latencias de uno a tres dias.

O dandole la vuelta con cuidado, tambien podriamos decir que usamos la solucion para dar una cota superior de r, dado R0 y el tiempo de latencia.

$$
{\ln R0 \over t_l} > r
$$

Autocompletar en python: Tries/Dagw vs suffix-arrays y familia

Esto es un problema que nos ha salido recientemente, tenemos que hacer sugerencias de autocompletar en las que no solo queremos los prefijos sino la posibilidad de completar cualquier match. Asi que el candidato obvio es alguna version avanzada del array de sufijos. Pues bien, no encuentro en Python 3 una libreria que sea lo suficientemente convincente.

He hecho unas cuantas benchmarks con un haystack de 18525 frases con longitud media de unos 22.3 caracteres. Ni muy grande ni muy pequeño, o tirando a pequeño, y por eso quizas es por lo que no hay forma de decidirse. He probado con 3 needles, una de un caracter, otra de 5 y otra de 11.

El tiempo a batir es el de una busqueda ingenua (NS), digamos [x for x in wordlist if subs in x], que tarda 7.6 milisegundos.

Estos han sido mis resultados

       1        5         11           
ASI    0.740   0.00428   0.00302
NS     7.35    7.62      7.63
dawg   3.63    0.00473   0.00157
datrie 7.20    0.00605   0.000947
reS    2.72    0.374     0.82
findS  2.68    0.542     0.404
pss     -      0.0539    0.0549
cst    4.15    0.00516   0.00315
STree  8.58    0.0252    0.0269

Los packages pysubstringsearch (pss) y suffix_trees (STree) se quedan claramente detras del package csuffixtree, como era de esperar. Pero tambien se quedan detras de un Array de Sufijos Ingenuo (ASI) que pone todas las frases en una lista ordenada y luego se limita a buscarla con bisect

[[KwList[binpos[n]] for n in range(bisect_left(binListData,kwsubs),bisect_left(binListData,kwsubs+"|||"))]

Otra idea ingenua es unir todo el texto y buscarlo con re o find (reS, findS). En este caso, a pesar de que no es muy grande el total, se nota la ausencia de un ordenado previo.

No veo pues otra forma de competir con el bisect que no sea usando una libreria subyacente en C o en algun lenguaje compilado. He probado los packages datrie y dawg, que en principio son para la prediccion de prefijos, pero podemos introducir a lo bruto las mismas claves del array de sufijos, esto es 391398 entradas. Son de hecho mas largas que en el array ya que para evitar repeticiones, lo que hago es poner rotaciones en vez de simplemente sufijos. A base de gastar mas memoria (más de la que se deberia para un Trie que aceptara repeticiones) por lo menos se consiguen buenas velocidades para needles altas; pero no tengo claro si merece la pena añadir una dependencia mas.