The following table presents data for the hemoglobin saturation dependin on \(O_2\) pressure:
mydata <- data.frame(
O2 = c(3.08, 4.61, 6.77, 10.15, 12.31, 15.38, 18.77, 22.77, 25.85, 30.15, 36.00, 45.23, 51.69, 61.85, 75.38, 87.08, 110.5),
Y = c(2.21, 3.59, 6.08, 10.50, 14.09, 19.34, 28.45, 40.33, 50.00, 60.50, 69.89, 80.11, 83.98, 88.95, 93.37, 95.86, 98.07)
)
knitr::kable(mydata, align = c("c", "c"))| O2 | Y |
|---|---|
| 3.08 | 2.21 |
| 4.61 | 3.59 |
| 6.77 | 6.08 |
| 10.15 | 10.50 |
| 12.31 | 14.09 |
| 15.38 | 19.34 |
| 18.77 | 28.45 |
| 22.77 | 40.33 |
| 25.85 | 50.00 |
| 30.15 | 60.50 |
| 36.00 | 69.89 |
| 45.23 | 80.11 |
| 51.69 | 83.98 |
| 61.85 | 88.95 |
| 75.38 | 93.37 |
| 87.08 | 95.86 |
| 110.50 | 98.07 |
Fit these data to a curve of the form
\[ \begin{align} Y = 100\cdot\frac{\sum_{j = 0}^4 j\alpha_j [O_2]^j}{4\sum_{j = 0}^4 \alpha_j [O_2]^j}. \end{align} \] Instruction: Suppose we have data points \(\{x_i, y_i\}\), \(i = 1, \dots, n\), that we wish to fit to some function \(y = f(x)\), and that the function \(f\) depends on parameters \(\{\alpha_i\}\), \(i = 1, \dots, m\). A fit of the data is achieved when the parameters are picked such that the function \[ \begin{align} F = \sum_{i = 1}^n (f(x_i) - y_i)^2 \end{align} \] is minimized. Fo find this fit, start with reasonable estimates for the parameters and then allow them to change dynamically (as a function of a time-like variable \(t\)) according to \[ \begin{align} \forall j:\quad \frac{d\alpha_j}{dt} = - \sum_{i = 1}^n \big(f(x_i)-y_i\big)\frac{\partial f(x_i)}{\partial \alpha_j}.\qquad\textrm{(1)} \end{align} \] With this choice, \[ \begin{align} \frac{dF}{dt} = \sum_{j = 1}^{m}\sum_{i = 1}^n 2\big(f(x_i)-y_i\big)\frac{\partial f(x_i)}{\partial \alpha_j}\frac{d\alpha_j}{dt}\leq 0, \end{align} \] so that \(F\) is a decreasing function of \(t\). A fit is found when numerical integration of \(\frac{d\vec\alpha}{dt}\) reaches a steady-state solution.
Implementation:
Define a function model(x, pars) that computes the model prediction \(Y\) from oxygen levels \([O_2]\), denoted by x and values \(\alpha_j\), denoted by pars. Compute derivatives \(\frac{\partial Y_i}{\partial \alpha_j}\) within this function and return the matrix of derivatives together with the model prediction vector, e.g. as an attribute or in a list. Advise: Write a small function to check whether the analytically computed derivates are correct.
Define a function func(times, states, parameters, mydata) for the dynamic system, eq. (1), to be used with ode(). Integrate the system and show the integration result on log-time scale. Check, if the integration time was sufficient to reach the steady-state.
Generate a plot showing the data points, the model prediction prior to optimization and after optimization.
What is the maximal degree \(j_{\max}\) needed to fit the data accurately? What is the measurement uncertainty of the data?
Episcopal churches have two towers. Freiburg is a diocesan town. Why does the cathedral have only one tower?