Mathematical tools for natural sciences

In this section we will learn to regress non-linear functions through data points. By non-liner we mean a function involving powers of x higher or lower than 1 and other functions like exponential and logarithmic functions.

There are two important non-linear functions used in many branches of sciences.

In the discussion below, we will learn to do regression using the method of least square fit to these two family of functions.

We need not struggle with a least square fit to this function. Our life is made simple by the fact that the log transformation of above relation can convert the exponential relation to a linear one!.

We have, \(~~\small{y = Ae^{Bx}~~}\), where A and B are constants.

Taking a natural logarithm on both sides,

\(\small{ log(y)~=~ log(A)~+~Bx }\)

Let \(\small{Y = log(y),~~~~b_0=log(A),~~~~and~~~b_1=B }\)

With this transformation, the exponential function \(\small{~~y = A e^{Bx} }\) is transformed to a linear function of the form,

convert \(\small{x_i, y_i}\) pairs of observations into, (x_i, Y_i=log(y_i)) pairs.

Then, using norml linear regression procedure we studied before, perform a least square fit to the straight line

Get back the original coefficients using,

\(\small{A = e^{b_0}~~~}\) and \(~~~\small{B = b_1}\)

Plot the fitted curve \(\small{x_i, Ae^{Bx_i}}\) on the same plot as aline with the data points \(\small{x_i,y_i }\)

The expected values of y is written as a polynomial of degree k as,

\(\small{y_{exp} = b_0 + b_1x + b_1x^2 + b_1x^3 + .... + b_kx^k }\)

We write down the difference between the expected and observed values summed over data points(called "residual"):

$$ \chi^2~=\sum_{i=1}^{n} [ y_i~-~(b_0+b_1x+b_2x^2+....+b_kx^k) ]^2 $$

We minimize the \(\small{\chi^2}\) values with respect to each one of the variables $b_0,b_1,b_2,....,b_k$ and equate them to zero.

\(\small{ \dfrac{\partial \chi^2}{\partial b_0}~=~ -2\displaystyle\sum_{i=1}^{n}[ y_i~-~(b_0+b_1x+b_2x^2+....+b_kx^k) ]~=~0 }\)

\(\small{ \dfrac{\partial \chi^2}{\partial b_1}~=~ -2\displaystyle\sum_{i=1}^{n}[ y_i~-~(b_0+b_1x+b_2x^2+....+b_kx^k) ] x~=~0 }\)

\(\small{ \dfrac{\partial \chi^2}{\partial b_2}~=~ -2\displaystyle\sum_{i=1}^{n}[ y_i~-~(b_0+b_1x+b_2x^2+....+b_kx^k) ]x^2~=~0 }\)

\(\small{ \dfrac{\partial \chi^2}{\partial b_3}~=~ -2\displaystyle\sum_{i=1}^{n}[ y_i~-~(b_0+b_1x+b_2x^2+....+b_kx^k) ]x^3~=~0 }\)

................................

.................................

\(\small{ \dfrac{\partial \chi^2}{\partial b_k}~=~ -2\displaystyle\sum_{i=1}^{n}[ y_i~-~(b_0+b_1x+b_2x^2+....+b_kx^k) ]x^k~=~0 }\)

This leads to a set of simultaneous equations given by,

\(\small{b_0n~+~ b_1\displaystyle\sum_{i=1}^{n}x_i~+~b_2\displaystyle\sum_{i=1}^{n}x_i^2~+~....~+b_k\displaystyle\sum_{i=1}^{n}x^k~=~\displaystyle\sum_{i=1}^{n} y_i }\)

\(\small{b_0\displaystyle\sum_{i=1}^{n}x_i~+~ b_1\displaystyle\sum_{i=1}^{n}x_i^2~+~b_2\displaystyle\sum_{i=1}^{n}x_i^3~+~....~+b_k\displaystyle\sum_{i=1}^{n}x_i^k~=~\displaystyle\sum_{i=1}^{n} y_i x_i }\)

........................................

.........................................

\(\small{b_0\displaystyle\sum_{i=1}^{n}x_i^k~+~ b_1\displaystyle\sum_{i=1}^{n}x_i^{k+1}~+~b_2\displaystyle\sum_{i=1}^{n}x_i^{k+2}~+~....~+b_k\displaystyle\sum_{i=1}^{n}x_i^{2k}~=~\displaystyle\sum_{i=1}^{n} x_i^ky_i }\)

The above set of k+1 simultaneous equations in variables \(\small{b_0, b_1,b_2,...,b_k }\) can be solved by matrix methods to get the coefficients \(\small{b_0, b_1,b_2,...,b_k }\).

Once the coefficients are obtained, the observed and fitted lines can be plotted on the same graph and $r^2$ value can be obtained, using methods similar to linear regression. We demonstrate how to do non-linear regression in R:

In the example R script below, we demonstrate lm() and nsl() functions of R library to fit the data.

# polynomial fit using lm ## define the two data sets x = c(0.5, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0) y = c(0.87, 7.12, 14.01, 24.37, 39.058, 58.62, 83.92) ## call lm() function with third degree polinomial of the form y = ax + bx^2 + cx^3 result = lm(y ~ x + I(x^2) + I(x^3)) # accessing coefficient of x print(result$coefficient[[1]]) # accessing coefficient of x^2 print(result$coefficient[[2]]) # accessing coefficient of x^3 print(result$coefficient[[3]]) #plot original points and the fitted curve plot(x,y,type="p", col="blue") yfitted = result$coefficient[[1]]*x + result$coefficient[[2]]*x^2 + result$coefficient[[3]]*x^3 lines(x,yfitted, type="l",col="red") # non-linear curve fitting with nls function # ---------------------------------------------------------------------------- # non-linear data is fitted with nls() function, by method of least squares. # we fit a power model. xdat = c(1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0) ydat = c(-0.09, 0.59, 0.42, 2.03, 3.43, 1.84, 5.30, 4.40, 6.67, 7.40, 7.34, 8.76, 10.25, 10.93, 13.78, 14.21, 17.82, 21.83, 23.04, 24.25, 27.48) df = data.frame(x = xdat, y = ydat) dat = nls(y ~ I(x^power), data=df) print(class(dat)) print(summary(dat)) print(summary(dat)$coefficients) power = round(summary(dat)$coefficients[1], 3) power.se = round(summary(dat)$coefficients[2], 3) s = seq(1, 3, 0.1) plot(ydat ~ xdat, main = "Fitted power model", sub = "Blue: fit; green: known") lines(s, s^3, lty = 2, col = "green") # We also fit a formula data = lm(formula = log(xdat)~ydat) summary(data)