Реализация метода наименьших квадратов в пакете Scilab

Пусть в результате эксперимента были получены некоторые данные (x_i, y_i) отображенные в виде таблицы. Требуется построить аналитическую зависимость

Y = f (x, a ₀ a ₁,..., a_k)

таким образом, чтобы сумма квадратов отклонений измеренных значений y_i от расчетных Y_i была наименьшей:

,

Задача сводится к определению коэффициентов a_i.

Для реализации этой задачи в Scilab предусмотрена функция

[a,S]=datafit(F,z,a0),

где

F – функция, с помощью которой можно вычислить разницу между измеренными и расчетными значениями: y - f (x, a ₀ a ₁,..., a_k);

z – матрица исходных данных (первая строка – массив x_i, вторая строка – y_i);

a0 – вектор начальных приближений (нулевой вектор из k элементов);

a – вектор коэффициентов (искомый вектор a ₀, a ₁,..., a_k);

S – сумма квадратов отклонений измеренных значений от расчетных.

Рассмотрим использование функции datafit на примере.

ПРИМЕР 4. Экспериментальные данные представлены в виде таблицы:

x

y

Выполнить регрессионный анализ. Методом наименьших квадратов подобрать зависимость вида y = c ₁ + c ₂ x + c ₃ x ² + c ₄ x ³.

Решение:

1. Ввод и графическое изображение исходных данных (рис. 18).

//Исходные данные

x=[1.32 1.40 1.50 1.62 1.70 1.80 1.90

2.00 2.11 2.20 2.32 2.40 2.51];

y=[3.30 3.50 3.85 4.25 4.50 4.85 5.40

6.00 6.60 7.30 9.20 10.20 13.50];

//Построение графика экспериментальных данных

plot2d(x,y,-4)

Рис. 18

2. Формирование матрицы исходных данных и вектора начальных приближений.

//Матрица исходных данных

z=[x;y];

3. Расчет коэффициентов линии регрессии y = a ₁ + a ₂ x, суммарной ошибки и коэффициента корреляции.

//Вектор начальных приближений

a=[0;0];

Разница между измеренными и расчетными значениями составляет:

y – (a ₁ + a ₂ x) = y – a ₁ – a ₂ x

В нашем случае x это первая строка матрицы z, y – вторая. Создадим функцию:

//Функция, задающая разницу между экспериментальными

//и теоретическими данными для линии регрессии.

function [zr]=R(a,z)

zr=z(2)-a(1)-a(2)*z(1)

endfunction

//Вычисление коэффициентов линии регрессии a

//и суммарной ошибки err

[a,err]=datafit(R,z,a)

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

//Результаты вычислений:

err =

15.526638

a =

- 7.3892172

7.2037863

Линия регрессии имеет вид: y = a ₁ + a ₂ x Þ y = -7.389 + 7.204 x.

Коэффициент корреляции вычисляется по формуле (8).

//Коэффициент корреляции

r=sum((x-mean(x)).*(y-mean(y)))/

sqrt(sum((x-mean(x)).^2)*sum((y-mean(y)).^2))

//Значение коэффициента корреляции

r =

0.9266121

4. Расчет коэффициентов линии y = c ₁ + c ₂ x + c ₃ x ² + c ₄ x ³, суммарной ошибки и индекса корреляции.

//Вектор начальных приближений

с=[0;0;0;0];

Разница между измеренными и расчетными значениями составляет:

y – (c ₁ + c ₂ x + c ₃ x ² + c ₄ x ³ ) = y – c ₁ – c ₂ x – c ₃ x ² – c ₄ x ³

В нашем случае x это первая строка матрицы z, y – вторая. Создадим функцию:

//Функция, задающая разницу между экспериментальными

//и теоретическими данными

function [zr]=G(c,z)

zr=z(2)-c(1)-c(2)*z(1)-c(3)*z(1)^2-c(4)*z(1)^3

endfunction

//Вычисление коэффициентов заданной линии с

//и суммарной ошибки err

[c,err]=datafit(G,z,c)

//Результаты вычислений:

err =

0.5287901

c =

- 51.576664

95.594671

- 55.695312

11.111453

Искомая линия имеет вид:

y = c ₁ + c ₂ x + c ₃ x ² + c ₄ x ³ Þ y = - 51.577+ 95.595 x + - 55.695 x ² + 11.111 x ³.

Расчет индекса корреляции выполняется по формуле (9).

//Индекс корреляции

R=sqrt(1-err/sum((y-mean(y))^2))

//Результаты вычислений:

R =

0.9975894

5. Графическое изображение поставленной задачи (рис. 19).

//Построение графиков экспериментальных данных и

//подобранных зависимостей в одной графической области

X=1.32:0.01:2.51;

YR=a(1)+a(2)*X;

YG=c(1)+c(2)*X+c(3)*X^2+c(4)*X^3;

plot(X,YR,X,YG)

Рис. 19

Дата добавления: 2015-12-08