Нелинейная регрессия

В своей работе химика-кинетика я часто встречался с задачей описать набор данных заданной нелинейной функцией. В R это делается достаточно просто.

1. Загружаем наши данные (предположим, что они содержатся в файле 1.csv):
> test=read.csv('~/1.csv',header=T)

2. На всякий случай проверяем, что получилось:
> test
X Y
1 0.5 0.00
2 1.0 0.01
3 2.0 0.02
4 3.0 0.05
5 4.0 0.13
6 5.0 0.18
7 6.0 0.24
8 7.0 0.29
9 8.0 0.35
10 9.0 0.37
11 10.0 0.41
12 11.0 0.44
13 12.0 0.47
14 15.0 0.57
15 20.0 0.71
16 25.0 0.80
17 30.0 0.86
18 35.0 0.90
19 43.0 0.95
20 50.0 0.97
21 60.0 0.99
22 70.0 0.99
23 80.0 1.00

Переменная X — время в минутах, Y — объём поглощённого кислорода в мл.

3. Попробуем построить график:
> attach(test)
> plot(Y~X)

4. Данные должны описываться кинетическим уравнением первого порядка Y = a*(1-exp(-k*X)), пробуем:
> test.mod=nls(Y~a*(1-exp(-k*X)),data=test,start=list(a=1,k=0.05))

Здесь nls() — команда для выполнения нелинейной регрессии, атрибутами которой являются:
а) модель, в нашем случае это Y~a*(1-exp(-k*X))
б) экспериментальные данные, у нас это test
в) стартовые значения параметров модели, которые подбираются исследователем эмпирически.

Если всё сделано правильно и параметры подобраны удачно, то в test.mod будет содержаться информация о полученной модели:
> test.mod
Nonlinear regression model
model: Y ~ a * (1 - exp(-k * X))
data: test
a k
1.05594 0.04967
residual sum-of-squares: 0.03621

Number of iterations to convergence: 4
Achieved convergence tolerance: 4.59e-06

5. Попробуем визуализировать нашу модель:
> test.predict=predict(test.mod,newdata=data.frame(X=0:100))

Здесь мы на промежутке от 0 минут до 100 вычисляем Y

Рисуем:
> plot(Y~X)
> lines(0:100,test.predict)

Ура! Всё получилось!