Context
Last week, I was discussing about how to use nls() for a specific
model with one of my colleague and I ended creating a piece of code to
show what I was talking about! Even though there are many posts exploring nls()
in more depth that I did (for instance this post on datascienceplus by Lionel Herzog),
I thought I could share these lines of command here!
Basically, we were talking about a model where the temperature ($T$) follows
a saturation curve starting from 10°C at t=0 (so T(0) = 10) and plateauing at \(T_{\inf}\).
$$T(t) = T_{\inf} - (T_{\inf} - T_0)\exp(-kt)$$
Goal and data
The goal here is to use nls() (Nonlinear Least Squares) to find \(k\) and \(T_{inf}\).
For the sack of clarity, I simulate the data, i.e. I use the saturation curve
with known parameter values, then I add some noise (here a white noise):
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
|
library(magrittr)
# Parameters
## known
T0 = 10
## the ones we are looking for
k = 0.2
Tinf = 20
# Simulate data
## time
seqt <- seq(0, 50, .25)
## create a data frame
simdata <- cbind(
seqt = seqt,
sim = Tinf - (Tinf- T0) * exp(-k * seqt) + .5 * rnorm(length(seqt))
) %>% as.data.frame
head(simdata)
#R> seqt sim
#R> 1 0.00 9.666048
#R> 2 0.25 12.027247
#R> 3 0.50 10.286901
#R> 4 0.75 11.167841
#R> 5 1.00 11.852154
#R> 6 1.25 12.653063
|
Use nls()
Now I call nls() to fit the data:
1
|
res <- nls(sim ~ Tinf - (Tinf - 10)*exp(-k*seqt), simdata, list(Tinf = 1, k = .1))
|
All the information needed are stored in res and display via the print method:
1
2
3
4
5
6
7
8
9
10
|
res
#R> Nonlinear regression model
#R> model: sim ~ Tinf - (Tinf - 10) * exp(-k * seqt)
#R> data: simdata
#R> Tinf k
#R> 19.9187 0.2094
#R> residual sum-of-squares: 52.4
#R>
#R> Number of iterations to convergence: 9
#R> Achieved convergence tolerance: 1.953e-06
|
Let’s draw a quick plot:
1
2
3
4
5
|
## get the coefficients values
cr <- coef(res)
fitC <- function(x) cr[1] - (cr[1] - 10)*exp(-cr[2]*x)
plot(simdata[,1], simdata[,2], xlab = "Time (min)", ylab = "Temperature (°C)")
curve(fitC, 0, 50, add = TRUE, col = "#2f85a4", lwd = 3)
|
That’s all folks!
Session info
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
|
sessionInfo()
#R> R version 4.2.2 (2022-10-31)
#R> Platform: x86_64-pc-linux-gnu (64-bit)
#R> Running under: Ubuntu 22.04.1 LTS
#R>
#R> Matrix products: default
#R> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#R> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
#R>
#R> locale:
#R> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#R> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#R> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#R> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#R>
#R> attached base packages:
#R> [1] stats graphics grDevices utils datasets methods base
#R>
#R> other attached packages:
#R> [1] magrittr_2.0.3 inSilecoRef_0.0.1.9000
#R>
#R> loaded via a namespace (and not attached):
#R> [1] tidyselect_1.2.0 xfun_0.35 bslib_0.4.1 vctrs_0.5.1
#R> [5] generics_0.1.3 miniUI_0.1.1.1 htmltools_0.5.4 yaml_2.3.6
#R> [9] utf8_1.2.2 rlang_1.0.6 later_1.3.0 pillar_1.8.1
#R> [13] jquerylib_0.1.4 httpcode_0.3.0 glue_1.6.2 DBI_1.1.3
#R> [17] lifecycle_1.0.3 plyr_1.8.8 stringr_1.5.0 blogdown_1.15
#R> [21] htmlwidgets_1.5.4 evaluate_0.18 knitr_1.41 fastmap_1.1.0
#R> [25] httpuv_1.6.6 curl_4.3.3 fansi_1.0.3 highr_0.9
#R> [29] Rcpp_1.0.9 xtable_1.8-4 promises_1.2.0.1 backports_1.4.1
#R> [33] DT_0.26 cachem_1.0.6 jsonlite_1.8.4 rcrossref_1.2.0
#R> [37] mime_0.12 digest_0.6.30 stringi_1.7.8 bookdown_0.30
#R> [41] dplyr_1.0.10 shiny_1.7.3 bibtex_0.5.0 cli_3.4.1
#R> [45] tools_4.2.2 sass_0.4.4 tibble_3.1.8 RefManageR_1.4.0
#R> [49] crul_1.3 pkgconfig_2.0.3 ellipsis_0.3.2 xml2_1.3.3
#R> [53] timechange_0.1.1 lubridate_1.9.0 assertthat_0.2.1 rmarkdown_2.18
#R> [57] httr_1.4.4 R6_2.5.1 compiler_4.2.2
|
