How to optimize the parameter values?

229 Views Asked by At

Let us assume that there is the number of infectious individuals, $I(t)$, where $t=1, 2, \ldots, T$ is time. We want to take into account the restrictive measures that affect the infection probability.

We introduced into a compartmental model, a scalar function of parameter for the restrictive measures: $$p(t): t \rightarrow [0, 1].$$ The step function $p(t)$ is a constant between two neighboring measures in time. Currently, the function's value defined by an expert in most cases after modeling:

\begin{equation*} p(t) \in \begin{cases} [0, 0.5), & \text{a restrictive measure limits the } I(t)\\ [0.5, 1], &\text{otherewise.} \end{cases} \end{equation*} Finally, we have obtained the model data $\hat I(t, p(t))$.

In the proposed model modification, we use four compartments with a cycle: Susceptible -> Exposed -> Infected -> Recovered -> Susceptible. Thus, we allow the recurrence of the disease and that's why we named the model SEIRS.

The figure shows the empirical data, $I(t)$, (red line) and model data, $\hat I(t, p(t))$, (blue dotted line), and the function, $p(t)$ (green line) where numbers indicate the function's values:

$$p(t)=\{0.29, 0.22, 0.17, 0.2, 0.24, 0.2, 0.17, 0.2, 0,25, 0.22, 0.17, 0.22\}.$$

One can see $0<p(t)\leq 0.29<0.5$ all times.

enter image description here

Question: it is required to optimize the function values, $p(t)$.

Edit after @Jean-ClaudeArbaut's comment

My attampt is:

  1. I have used the approximation accuracy as the optimality criterion: $$|I(t)-\hat I(t, p(t))|$$ then a target function can be $$L(t, p(t))=\sqrt{\frac{1}{T} \sum_{t=1}^T(I(t)-\hat I(t, p(t)))^2} \rightarrow \underset{p(t)}\min \\ \text{ s. t. } p(t) \in [0, 1]\\ t = 1, 2, \ldots T.$$
  1. I have took the $p_1=0.29$, and setup the lower and upper bounds, then tried minimizing $L$ over $p_1$ for first time window where $t=1,2, \ldots, 50$ with parameter fitting using Levenberg-Marquart algorithm:
    library(minpack.lm) 
    
    I <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 26, 32, 45, 
    64, 88, 122, 136, 154, 165, 176, 192, 214, 232, 254, 
    275, 301, 338, 360, 395, 400, 413, 440, 453, 467, 474, 
    498, 515, 538, 550, 570, 598, 614, 611, 618, 619, 614, 
    614, 616, 600))
    
    I_hat <-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 28, 37, 52, 
    74, 98, 133, 150, 163, 178, 201, 233, 261, 282, 310, 
    334, 359, 382, 413, 439, 439, 460, 483, 494, 506, 495, 
    509, 527, 541, 553, 567, 586, 609, 609, 602, 586, 569, 
    560, 572, 572)
    
    ## residual function 
    residual <- function(I_hat, observed) I_hat - observed 
    
    # parameter fitting using Levenberg-Marquart algorithm
     fitval = nls.lm(par = 0.29, lower=0, upper=1, fn = residual, 
                     observed = I, control = nls.lm.control(nprint=1))

The output is:

It.    0, RSS = 6.85828e+06, Par. =       0.29
It.    1, RSS = 6.8379e+06, Par. =          1
It.    2, RSS = 4.73522e+06, Par. =          1

One can see that $p_1$ obtains the max value on the 1st step.