A Fit of Moscow Covid-19 Data

  In the last blog it was shown that knowing I_max, the peak value for the number of Infected individuals, helped reduce the number of parameters needed to estimate the values of the other parameters for a curve fit. In the case of a Covid-19 fit the peak value is not always immediately known. At the peak we know that dI/dt=(rS-a)I=0 which implies that rS-a=0 but we don't know what the values of S are and only have values of I to work with. However we can evaluate d(lnI)/dt=rS-a and look for a zero of this function.

In the above plot of Moscow I data we see that this function is about to cross the zero horizontal axis if it hasn't done so already. To get an estimated set of parameters for the fit the value of  Imax was arbitrarily assumed to be 11000. Then a search of the remaining parameters was done to minimize the value of the rms error for the fit.

New data will allow us to improve on this fit and get a better estimate of the maximum value of I.

Reducing the Set of Parameters for the Fit

  One needs a set of 5 parameters such as N, I₀, S₀, r, and a to do the influenza curve fit. Since N=763 is given and I₀ can be taken to be roughly the value for the 1st data point we are left with 3 unknowns for a rough fit. If the data includes the peak of the Infected values we can use it and a formula for I_max in terms of the other unknowns to get an estimate for S₀. We are then left with just 2 unknowns needed to make a rough estimate of the parameters which we will take to be r and ρ=a/r.

Given an estimate for I_max one can compute an table of values for S₀(ρ), then fit a quadratic curve to the data from the table and use the quadratic to interpolate the data for an assumed value of ρ to estimate S₀.

Once a rough estimate is made one can adjust the set of 4 parameters to give a minimum value for the rms error.

The value of "δI rms err" is the rms error for the difference between I(fit) and I(formula). I added a switch to the spreadsheet to allow a comparison with a 1st order numerical integration.

 The 2nd order calculation is definitely better.

A Sample Influenza Curve Fit

  A curve fit can be somewhat challenging if you don't have a formula to actually fit to your data but it's not impossible if you have a set of differential equations that model the process. I've worked out a procedure that gives fairly good results for the SIR epidemic model and applied it to the 1978 English boys school influenza epidemic found in Murray, Mathematical Biology. The set of equations are shown in the following image along with an integrated expression for the number of infected, I, as a function of the unknown number of susceptibles, S. The 3rd compartment for the SIR model is the number of removed individuals, R, but we will not need that.

Given assumed initial values for I and S one can compute successive values using a simple numerical integration procedure. A first order calculation appears to be sufficient if the step size is small enough. I tried a second order term to help get past the peak but it didn't make much difference in the results.

The data used for the fit was obtained with the aid of verinier calipers to measure the height of the data points in the figure found in Murray's book.

The small step size chosen required about 600 iterations for the numerical integration to cover the range of the data. The fits were evaluated by comparing the rms errors of the fits and values for the set of parameters of the model were chosen by trial and error to obtain the smallest rms error.

Here's a plot showing a comparison of the fit with the data along with a comparison of the values of I(S) for both the numerical integration and the formula which was used as a check on the accuracy of the numerical integration.

Some other plots are also useful.

Supplemental (May 15): The 2nd plot from the end should be labeled I(S).

Using Google Sheets to Check the Fit

    Since I was limited to the use my mobile phone I installed Google Sheets to do some simple calculations. Below is a check on the fit for the cooling temperature data and constants given in the last two blogs. Just below constants are three lines that show the formulas for the indicated cells.
This fit shows a slight improvement over the original fit due to a change in the interval used to fit the thermistor temperature curve as indicated by the cell containing the rms error. The rounded values for the temperatures may have also affected the fit.

Missing Cooling Experiment Rate Constants

  I inadvertantly failed to include the two rate constants, λ1 and λ2, for the fit of the cooling data in the last blog. One can now use the constants given here to check the fit.

The fit was rather difficult to do since the formula for the fit function has a nonlinear set of undetermined constants so one cannot use linear least squares. The solution required a search for the values of the rate constants along with a linear fit of the constants A, B and C that would minimize the rms error of the fit. In addition, the problem was ill-conditioned.

A More Precise Cooling Experiment

  Two weeks ago I redid the cooling experiment paying more attention to the precision measurement of the temperature through the use of a thermistor. The experimental setup is shown below with the thermistor inserted in a tall glass of hot tap water.

The thermistor was insulated with heat shrink tubing and held in place by a foam board card. The thermistor resistance was measured with a digital multimeter. The thermistor was obtained from Radio Shack a few years ago and contained calibration data on the back of the card.

A cubic fit of the characteristic curve for the temperature vs the natural log of the relative resistance allowed the resistance measurements collected to be converted to temperatures in °C. The electrical properties of a thermistor can be found in Millman & Halkias, Electronic devices and circuits, 1967, McGraw-Hill, p.132.

The resulting measurements were very consistent.

Using an electrical analogy for the cooling mechanism one can derive a formula involving exponential functions that one can fit to the data. The fit is quite good and also compares well with the rate of internal cooling which is proportional to the negative rate of change in temperature.

  The data for the last plot used a numerical formula to estimate the rate of change in temperature.

The results for the experiment were quite good but that doesn't necessarily mean they were accurate since we couldn't check thermistor calibration against a standard gas thermometer or, equivalently, the thermodynamic temperature scale. The ambient temperature of the room remained within a few tenths of 70.1 °F throughout the entire time the data was collected and there may have been a some slight warming towards the end of the data. Although I recorded the data manually the entire experiment could be done with a computer recording the data using a USB port if the appropriate drivers are installed on the computer.

One wonders how accurately the temperatures are measured using weather temperature recording programs?

An Outline of a Derivation of an Algebraic Expression for cos(1°)

  Yesterday I posted a tweet containing a closed-form algebraic expression for the cosine of 1°. One might wonder how one can arrive at this formula. The answer is that one has to use the formulas of trigonometry for combining angles, or the formulas for the product of complex numbers or alternatively formulas for successive rotations. We'll start by giving a few formulas for rotations. The following matrix rotates a point on the x,y-plane an angle counter clockwise through an angle θ.

A number of useful formulas are the following in which a product is equated to the identity matrix which corresponds to a complete rotation of 360 degrees or another matrix for some other known rotation.

Notice that one can use the cube of a rotation matrix to get a cubic formula in reduced form involving only the cosine of the unknown angle and that of the known angle. One can use the Cardano formula to solve for κ. The last formula is for the difference between two known angles. These formulas can be used to find the sines and cosines of the following angles.

The solution of the first four matrix power formulas is simplified by using the components equal to zero and eliminating unneeded factors. One has to be careful about the κ in the fourth power equation since the solution is κ=0 for θ=90°. A check shows that the formulas used in the tweet work as desired.