Метод наименьших квадратов (МНК) позволяет по экспериментальным данным подобрать такую аналитическую функцию, которая проходит настолько близко к экспериментальным точкам, насколько это возможно.
В общем случае задачу можно сформулировать следующим образом. Пусть в результате эксперимента были получены некая экспериментальная зависимость $$y(x)$$, представленная в таблице 11.1. Необходимо построить аналитическую зависимость $$f(x,a_{1},a_{2},\dots, a_{k})$$, наиболее точно описывающую результаты эксперимента. Для построения параметров функции $$f(x,a_{1},a_{2},\dots, a_{k})$$ будем использовать метод наименьших квадратов. Идея метода наименьших квадратов заключается в том, что функцию $$f(x,a_{1},a_{2},\dots, a_{k})$$ необходимо подобрать таким образом, чтобы сумма квадратов отклонений измеренных значений $$Y_i=f(x,a_{1},a_{2},\dots, a_{k})$$ была бы наименьшей (см. рис. 11.1):
| x | $$x_1$$ | $$x_2$$ | $$x_3$$ | ... | $$x_n-1$$ | $$x_n$$ |
| y | $$y_1$$ | $$y_2$$ | $$y_3$$ | ... | $$y_n-1$$ | $$y_n$$ |
Задача состоит из двух этапов:
Математически задача подбора коэффициентов зависимости сводится к определению коэффициентов $$a_i$$ из условия (11.1). В Octave её можно решать несколькими способами:
Вспомним некоторые сведения из высшей математики, необходимые для решения задачи подбора зависимости методом наименьших квадратов.
Достаточным условием минимума функции $$S(a_{1},a_{2},\dots, a_{k})$$ является равенство нулю всех её частных производных. Поэтому задача поиска минимума функции (11.1) эквивалентна решению системы
(рис 11.1) Геометрическая интерпретация МНК
алгебраических уравнений:
$$\left\{ \begin{matrix} \frac{\partial S}{\partial a_1}=0\\ \frac{\partial S}{\partial a_2}=0\\ \dots\\ \frac{\partial S}{\partial a_k}=0 \end{matrix} \right.$$Если параметры $$a_i$$ входят в зависимость $$Y=f(x,a_{1},a_{2},\dots, a_{k})$$ линейно, то получим систему (11.3) из k линейных уравнений с $$k$$ неизвестными.
$$\left\{ \begin{matrix} \displaystyle \sum_{i=1}^n2\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)\frac{\partial f}{\partial a_1}=0\\ \displaystyle \sum_{i=1}^n2\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)\frac{\partial f}{\partial a_2}=0\\ \hdotsfor{1}\\ \displaystyle \sum_{i=1}^n2\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)\frac{\partial f}{\partial a_k}=0 \end{matrix} \right.$$Составим систему (11.3) для наиболее часто используемых функций.
Для подбора параметров линейной функции $$Y=a_1+a_2x$$, составим функцию (11.1) для линейной зависимости:
$$S(a_1,a_2)=\sum_{i=1}^n\left(y_i-a_1-a_2x_i\right)^2\to \min.$$Продифференцировав функцию $$S$$ по $$a_1$$ и $$a_2$$, получим систему уравнений:
$$\left\{ \begin{matrix} \displaystyle 2\sum_{i=1}^n\left(y_i-a_1-a_2x_i\right)(-1)=0\\ \displaystyle 2\sum_{i=1}^n\left(y_i-a_1-a_2x_i\right)(-x_i)=0 \end{matrix} \right. \Rightarrow \left\{ \begin{matrix} \displaystyle a_1n+a_2\sum_{i=1}^nx_i=\sum_{i=1}^ny_i\\ \displaystyle a_1\sum_{i=1}^nx_i+a_2\sum_{i=1}^nx_i^2=\sum_{i=1}^ny_ix_i \end{matrix} \right.$$решив которую, определим коэффициенты функции $$Y=a_1+a_2x$$:
$$\left\{ \begin{matrix} \displaystyle a_1=\displaystyle \frac{\displaystyle\sum_{i=1}^ny_i}n-a_2\displaystyle \frac{\displaystyle \sum_{i=1}^nx_i}n\\ \displaystyle a_2=\displaystyle \frac{n\displaystyle\sum_{i=1}^ny_ix_i-\sum_{i=1}^ny_i\sum_{i=1}^nx_i}{n\displaystyle\sum_{i=1}^nx_i^2-\left(\sum_{i=1}^nx_i\right)^2} \end{matrix} \right.$$Для определения параметров зависимости $$Y=a_1+a_2x+a_3x^2$$ составим функцию $$S(a_1,a_2,a_3)$$(11.1):
$$S(a_1,a_2,a_3)=\sum_{i=1}^n\left(y_i-a_1-a_2x_i-a_3x_i^2\right)^2\to \min$$После дифференцирования S по $$a_i$$,$$a_2$$ и $$a_3$$ получим систему линейных алгебраических уравнений:
$$\left\{ \begin{matrix} \displaystyle a_{1}n+a_{2}\sum_{i=1}^{n}x_{i}+a_{3}\sum_{i=1}^{n}x_{i}^{2}=\sum_{i=1}^{n}y_{i}\\ \displaystyle a_{1}\sum_{i=1}^{n}x_{i}+a_{2}\sum_{i=1}^{n}x_{i}^{2}+a_{3}\sum_{i=1}^{n}x_{i}^{3}=\sum_{i=1}^{n}y_{i}x_{i}\\ \displaystyle a_{1}\sum_{i=1}^{n}x_{i}^{2}+a_{2}\sum_{i=1}^{n}x_{i}^{3}+a_{3}\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}y_{i}x_{i}^{2} \end{matrix} \right.$$Решив систему (11.8), найдём значения параметров a$$a_i$$,$$a_2$$ и $$a_3$$.
Аналогично определим параметры многочлена третьей степени:$$S(a_1,a_2,a_3,a_4)$$. Составим функцию $$S(a_1,a_2,a_3,a_4)$$:
$$S(a_1,a_2,a_3,a_4)=\sum_{i=1}^n\left(y_i-a_1-a_2x_i-a_3x_i^2-a_4x_i^3\right)^2\to\min$$После дифференцирования $$S$$ по $$a_i$$,$$a_2$$,$$a_3$$ и $$a_4$$, система линейных алгебраических уравнений для вычисления параметров a $$a_1,a_2,a_3,a_4$$ примет вид:
$$\left\{ \begin{matrix} \displaystyle a_{1}n+a_{2}\sum_{i=1}^{n}x_{i}+a_{3}\sum_{i=1}^{n}x_{i}^{2}+a_{4}\sum_{i=1}^{n}x_{i}^{3}=\sum_{i=1}^{n}y_{i}\\ \displaystyle a_{1}\sum _{i=1}^{n}x_{i}+a_{2}\sum _{i=1}^{n}x_{i}^{2}+a_{3}\sum_{i=1}^{n}x_{i}^{3}+a_{4}\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}y_{i}x_{i}\\ \displaystyle a_{1}\sum _{i=1}^{n}x_{i}^{2}+a_{2}\sum_{i=1}^{n}x_{i}^{3}+a_{3}\sum _{i=1}^{n}x_{i}^{4}+a_{4}\sum_{i=1}^{n}x_{i}^{5}=\sum_{i=1}^{n}y_{i}x_{i}^{2}\\ \displaystyle a_{1}\sum_{i=1}^{n}x_{i}^{3}+a_{2}\sum_{i=1}^{n}x_{i}^{4}+a_{3}\sum_{i=1}^{n}x_{i}^{5}+a_{4}\sum_{i=1}^{n}x_{i}^{6}=\sum_{i=1}^{n}y_{i}x_{i}^{3} \end{matrix} \right.$$Решив систему (11.10), найдём коэффициенты $$a_i$$,$$a_2$$,$$a_3$$ и $$a_4$$.
В общем случае система уравнений для вычисления параметров $$a_i$$ многочлена k-й степени $$Y=\sum_{i=1}^{k+1}a_{i}x^{i-1}$$ имеет вид:
$$\left\{ \begin{matrix} \displaystyle a_1n+a_2\sum_{i=1}^nx_i+a_3\sum_{i=1}^nx_i^2+\ldots +a_{k+1}\sum_{i=1}^nx_i^k=\sum_{i=1}^ny_i\\ \displaystyle a_1\sum_{i=1}^nx_i+a_2\sum_{i=1}^nx_i^2+a_3\sum_{i=1}^nx_i^3+\ldots +a_{k+1}\sum_{i=1}^nx_i^{k+1}=\sum_{i=1}^ny_ix_i\\ \hdotsfor{1}\\ \displaystyle a_1\sum_{i=1}^nx_i^{k+1}+a_2\sum_{i=1}^nx_i^{k+2}+a_3\sum_{i=1}^nx_i^{k+3}+\ldots +a_{k+1}\sum_{i=1}^nx_i^{2k}=\sum_{i=1}^ny_ix_i^k \end{matrix} \right.$$В матричном виде систему (11.11) можно записать
$$Ca = g,$$Элементы матрицы C и вектора g рассчитываются по формулам
$$C_{i,j}=\sum_{k=1}^nx_k^{i+j-2},\quad i=1,\dots,k+1,\ j=1,\dots,k+1$$ $$g_i=\sum_{k=1}^ny_kx_k^{i-1},\quad i=1,\dots,k+1$$Решив систему (11.12), определим параметры зависимости $$Y=a_1+a_2x+a_3x^2+\dots+a_{k+1}x^k$$.
$$Y=ax^be^{cx}$$
Параметры $$b$$ и $$c$$ входят в зависимость $$Y=ax^be^{cx}$$ нелинейным образом. Чтобы избавиться от нелинейности предварительно
Составим функцию $$S(A, b, c)$$ по формуле (11.1):
$$S(A,b,c)=\sum_{i=1}^n\left(Y1_i-A-b\ln x_i-cx_i\right)^2\to\min$$После дифференцирования получим систему трёх линейных алгебраических уравнений для определения коэффициентов $$A, b, c$$.
$$\left\{ \begin{matrix} \displaystyle nA+b\sum_{i=1}^n\ln x_i+c\sum_{i=1}^nx_i=\sum_{i=1}^nY1_i\\ \displaystyle A\sum_{i=1}^n\ln x_i+b\sum_{i=1}^n(\ln x_i)^2+b\sum_{i=1}^nx_i\ln x_i=\sum_{i=1}^nY1_i\ln x_i\\ \displaystyle A\sum_{i=1}^nx_i+b\sum_{i=1}^nx_i\ln x_i+c\sum_{i=1}^nx_i^2=\sum_{i=1}^nY1_ix_i \end{matrix} \right.$$После решения системы (11.16) необходимо вычислить значение коэффициента a по формуле $$a=e^A$$.
Для вычисления параметров функции $$Y = ax^b$$ необходимо предварительно её прологарифмировать $$\ln Y=\ln ax^b=\ln a+b\ln x$$. После чего замена $$Z = lnY, X = lnx, A = lna$$ приводит заданную функцию к линейному виду $$Z = bX + A$$, где коэффициенты $$A$$ и $$b$$ вычисляются по формулам (11.6) и, соответственно, $$a=e^A$$.
Аналогично можно подобрать параметры функции вида $$Y = ae^{bx}$$. Прологарифмируем заданную функцию $$\ln y=\ln a+bx\ln e,\ln y=\ln a+bx$$. Проведём замену $$Y = lny, A = lna$$ и получим линейную зависимость $$Y = bx + A$$. По формулам (11.6) найдём $$A$$ и $$b$$, а затем вычислим $$a=e^A$$.
Рассмотрим ещё ряд зависимостей, которые сводятся к линейной.
Для подбора параметров функции $$Y=\frac{1}{ax+b}$$ сделаем замену $$Y=\frac{1}{Y}$$. В результате получим линейную зависимость Z = ax + b. Функция $$Y=\frac{x}{ax+b}$$ заменами $$Z=\frac{1}{Y}$$, $$X=\frac{1}{x}$$ сводится к линейной Z = a + bX. Для определения коэффициентов функциональной зависимости $$Y=\frac{1}{ae^{-x}+b}$$ необходимо сделать следующие замены $$Z=\frac{1}{Y},X=e^{-x}$$. В результате также получим линейную функцию Z = aX + b.
Аналогичными приёмами (логарифмированием, заменами и т. п.) можно многие подбираемые зависимости преобразовать к такому виду, что получаемая при решении задачи оптимизации система (11.2) была системой линейных алгебраических уравнений. При использовании Octave можно напрямую решать задачу подбора параметров, как задачу оптимизации (11.1) с использованием функции sqp.
После нахождения параметров зависимости $$f(x,a_{1},a_{2},\dots, a_{k})$$ возникает вопрос насколько адекватно описывает подобранная зависимость экспериментальные данные. Чем ближе величина
$$S=\sum_{i=1}^n\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)^2$$называемая суммарной квадратичной ошибкой, к нулю, тем точнее подобранная кривая описывает экспериментальные данные.
Линия, описываемая уравнением вида $$y=a_{1}+a_{2}x$$, называется линией регрессии $$y$$ на $$x$$, параметры $$a_1$$ и $$a_2$$ называются коэффициентами регрессии и определяются формулами (11.6).
Чем меньше величина $$S=\sum _{i=1}^{n}(y_{i}-a_{1}-a_{2}x_{i})^{2}$$, тем более обоснованно предположение, что экспериментальные данные описываются линейной функцией. Существует показатель, характеризующий тесноту линейной связи между $$x$$ и $$y$$, который называется коэффициентом корреляции и рассчитывается по формуле:
$$r=\frac{\displaystyle \sum_{i=1}^n\left(x_i-M_x\right)\left(y_i-M_y\right)} {\sqrt{\displaystyle\sum_{i=1}^n\left(x_i-M_x\right)^2\sum_{i=1}^n\left(y_i-M_y\right)^2}}, \ M_{x}=\frac{\displaystyle\sum_{i=1}^nx_i}{n},\ M_{y}=\frac{\displaystyle\sum_{i=1}^ny_i}{n}$$Значение коэффициента корреляции удовлетворяет соотношению $$-1\le r\le 1$$.
Чем меньше отличается абсолютная величина $$r$$ от единицы, тем ближе к линии регрессии располагаются экспериментальные точки. Если $$|r| = 1$$, то все экспериментальные точки находятся на линии регрессии. Если коэффициент корреляции близок к нулю, то это означает, что между $$x$$ и $$y$$ не существует линейной связи, но между ними может существовать зависимость, отличная от линейной.
Для того, чтобы проверить, значимо ли отличается от нуля коэффициент корреляции, можно использовать критерий Стьюдента. Вычисленное значение критерия определяется по формуле:
$$t=r\sqrt{\frac{n-2}{1-r^2}}$$Рассчитанное по формуле (11.19) значение $$t$$ сравнивается со значением, взятым из таблицы распределения Стьюдента (см. табл. 11.2) в соответствии с уровнем значимости $$p$$ (стандартное значение $$p = 0.95$$) и числом степеней свободы $$k = n - 2$$. Если полученная по формуле (11.19) величина $$t$$ больше табличного значения, то коэффициент корреляции значимо отличен от нуля.
| k\p | 0,99 | 0,98 | 0,95 | 0,90 | 0,80 | 0,70 | 0,60 |
| 1 | 63,657 | 31,821 | 12,706 | 6,314 | 3,078 | 1,963 | 1,376 |
| 2 | 9,925 | 6,965 | 4,303 | 2,920 | 1,886 | 1,386 | 1,061 |
| 3 | 5,841 | 4,541 | 3,182 | 2,353 | 1,638 | 1,250 | 0,978 |
| 4 | 4,604 | 3,747 | 2,776 | 2,132 | 1,533 | 1,190 | 0,941 |
| 5 | 4,032 | 3,365 | 2,571 | 2,05 | 1,476 | 1,156 | 0,920 |
| 6 | 3,707 | 3,141 | 2,447 | 1,943 | 1,440 | 1,134 | 0,906 |
| 7 | 3,499 | 2,998 | 2,365 | 1,895 | 1,415 | 1,119 | 0,896 |
| 8 | 3,355 | 2,896 | 2,896 | 1,860 | 1,387 | 1,108 | 0,889 |
| 9 | 3,250 | 2,821 | 2,261 | 1,833 | 1,383 | 1,100 | 0,883 |
| 10 | 3,169 | 2,764 | 2,228 | 1,812 | 1,372 | 1,093 | 0,879 |
| 11 | 3,106 | 2,718 | 2,201 | 1,796 | 1,363 | 1,088 | 0,876 |
| 12 | 3,055 | 2,681 | 2,179 | 1,782 | 1,356 | 1,083 | 0,873 |
| 13 | 3,012 | 2,681 | 2,179 | 1,782 | 1,356 | 1,083 | 0,873 |
| 14 | 2,977 | 2,624 | 2,145 | 1,761 | 1,345 | 1,076 | 0,868 |
| 15 | 2,947 | 2,602 | 2,131 | 1,753 | 1,341 | 1,074 | 0,866 |
| 16 | 2,921 | 2,583 | 2,120 | 1,746 | 1,337 | 1,071 | 0,865 |
| 17 | 2,898 | 2,567 | 2,110 | 1,740 | 1,333 | 1,069 | 0,863 |
| 18 | 2,878 | 2,552 | 2,101 | 1,734 | 1,330 | 1,067 | 0,862 |
| 19 | 2,861 | 2,539 | 2,093 | 1,729 | 1,328 | 1,066 | 0,861 |
| 20 | 2,845 | 2,528 | 2,086 | 1,725 | 1,325 | 1,064 | 0,860 |
| 21 | 2,831 | 2,518 | 2,080 | 1,721 | 1,323 | 1,063 | 0,859 |
| 22 | 2,819 | 2,508 | 2,074 | 1,717 | 1,321 | 1,061 | 0,858 |
| 23 | 2,807 | 2,500 | 2,069 | 1,714 | 1,319 | 1,060 | 0,858 |
| 24 | 2,797 | 2,492 | 2,064 | 1,711 | 1,318 | 1,059 | 0,857 |
| 25 | 2,779 | 2,485 | 2,060 | 1,708 | 1,316 | 1,058 | 0,856 |
| 26 | 2,771 | 2,479 | 2,056 | 1,706 | 1,315 | 1,058 | 0,856 |
| 27 | 2,763 | 2,473 | 2,052 | 1,703 | 1,314 | 1,057 | 0,855 |
| 28 | 2,756 | 2,467 | 2,048 | 1,701 | 1,313 | 1,056 | 0,855 |
| 29 | 2,750 | 2,462 | 2,045 | 1,699 | 1,311 | 1,055 | 0,854 |
| 30 | 2,704 | 2,457 | 2,042 | 1,697 | 1,310 | 1,055 | 0,854 |
| 40 | 2,660 | 2,423 | 2,021 | 1,684 | 1,303 | 1,050 | 0,851 |
| 60 | 2,612 | 2,390 | 2,000 | 1,671 | 1,296 | 1,046 | 0,848 |
| 120 | 2,617 | 2,358 | 1,980 | 1,980 | 1,289 | 1,041 | 0,845 |
| $$\infty$$ | 2,576 | 2,326 | 1,960 | 1,645 | 1,282 | 1,036 | 0,842 |
Коэффициент корреляции $$r$$ применяется только в тех случаях, когда между данными существует прямолинейная связь. Если же связь нелинейная, то для выявления тесноты связи между переменными $$y$$ и $$x$$ в случае нелинейной зависимости пользуются индекс корреляции. Он показывает тесноту связи между фактором $$x$$ и зависимой переменной $$y$$ и рассчитывается по формуле:
$$R=\sqrt{1-\frac{\displaystyle \sum _{i=1}^{n}\left(y_{i}-Y_{i}\right)^{2}}{\sum_{i=1}^{n}\left(y_{i}-M_{y}\right)^{2}}}$$где $$y$$ — экспериментальные значения, $$Y $$— теоретические значения (рассчитанные по подобранной методом наименьших квадратов формуле), $$M_y$$ — среднее значение $$y$$.
Индекс корреляции лежит в пределах от 0 до 1. При наличии функциональной зависимости индекс корреляции близок к 1. При отсутствии связи $$R$$ практически равен нулю. Если коэффициент корреляции r является мерой тесноты связи только для линейной формы связи, то индекс корреляции $$R$$ — как для линейной, так и для нелинейной. При прямолинейной связи коэффициент корреляции по своей абсолютной величине равен индексу корреляции: $$|r| = R$$.
| t | $$0^o$$ | $$4^o$$ | $$10^o$$ | $$15^o$$ | $$21^o$$ | $$29^o$$ | $$36^o$$ | $$51^o$$ | $$68^o$$ |
| P | -66,7 | -71,0 | 76,3 | 80,6 | 85,7 | 92,9 | 99,4 | 113,6 | 125,1 |
Для решения задач подбора аналитических зависимостей по экспериментальным данным можно использовать следующие функции Octave:
Пример 11.1. В ->Основах химии-> Д.И. Менделеева приводятся данные о растворимости азотнокислого натрия $$NaNO_3$$ в зависимости от температуры воды. В 100 частях воды (табл. 11.3) растворяется следующее число условных частей $$NaNO_3$$ при соответствующих температурах. Требуется определить растворимость азотнокислого натрия при температуре $$t = 32°C$$ в случае линейной зависимости и найти коэффициент корреляции.
Решение задачи 11.1 с комментариями приведено в листинге 11.1.
(рис 11.2) Иллюстрация к примеру 11.1
% Ввод экспериментальных данных X=[0 4 10 15 21 29 36 51 68]; Y=[66.7 71.0 76.3 80.6 85.7 92.9 99. 4 113.6 125.1]; % Вычисление вектора коэффициентов полинома y=a1*x+a2 [ a]= polyfit (X, Y, 1) % Вычисление значения полинома y=a1*x+a2 в точке t=32 t =32; yt=a(1)*t+a(2) % Построение графика полинома y=a1*x+a2, экспериментальных точек % и значения в заданной точке в одной графической области x =0:68; y=a(1)*x+a(2); plot (X, Y, ’ok’, x, y, ’-k’, t, yt, ’*k’) grid on; % Вычисление коэффициента корреляции k =cor(x, y) % Результаты вычислений a = 0.87064 67.50779 yt = 95.368 k = 1
На рис. 11.2 приведено графическое решение этой задачи, изображены экспериментальные точки, линия регрессии $$y=a_1x+a_2$$, на котором отмечена точка $$t = 32$$.
Пример 11.2. В результате эксперимента получена табличная зависимость $$y(x)$$ (см. табл. 11.4). Подобрать аналитическую зависимость $$Y=ax^be^{cx}$$ методом наименьших квадратов. Вычислить ожидаемое значение в точках 2, 3, 4. Вычислить индекс корреляции. Решение задачи подбора параметров функции $$f(x)=ax^be^{cx}$$ в Octave возможно двумя способами:
| x | 1 | 1,4 | 1,8 | 2,2 | 2,6 | 3 | 3,4 | 3,8 | 4,2 | 4,6 | 5 | 5,4 | 5,8 |
| y | 0,7 | 0,75 | 0,67 | 0,62 | 0,51 | 0,45 | 0,4 | 0,32 | 0,28 | 0,25 | 0,22 | 0,16 | 0,1 |
Рассмотрим последовательно оба варианта решения задачи.
Способ 1.
Функция (11.15) реализована в Octave с помощью функции $$f\_mnk$$. Полный текст программы решения задачи способом 1 с комментариями приведён на листинге 11.2. Вместо коэффициентов $$A, b, c$$ из формул (11.15)–(11.16) в программе на Octave используется массив $$c$$.
function s=f_mnk( c ) % Переменные x,y являются глобальными, используются в нескольких функциях global x; global y; s =0; for i =1: length(x) s=s +(log(y(i))-c(1)-c(2)*log(x(i))*c(3)*x(i))^2; end end % _________________________________________________ global x; global y; % Задание начального значения вектора c, при неправильном его % определении, экстремум может быть найден неправильно. c = [2; 1; 3]; % Определение координат экспериментальных точек x=[1 1.4 1.8 2.2 2.6 3 3.4 3.8 4.2 4.6 5 5.4 5.8]; y=[0.7 0.75 0.67 0.62 0.51 0.45 0.4 0.32 0.28 0.25 0.22 0.16 0.1]; % Решение задачи оптимизации функции 11.15 с помощью sqp. c=sqp(c,@f_mnk) % Вычисление суммарной квадратичной ошибки для подобранной % зависимости и вывод её на экран. sum1=f_mnk( c ) % Формирование точек для построения графика подобранной кривой. x1 = 1:0.1:6; y1= exp(c(1)).*x1.^ c(2).*exp(c(3).*x1); % Вычисление значений на подобранной кривой в заданных точках. yr=exp(c(1)).*x.^ c(2).*exp(c(3).*x ); % Вычисление ожидаемого значения подобранной функции в точках x=[2,3,4] x2 =[2 3 4 ] y2= exp(c(1)).*x2.^c(2).*exp(c(3).*x2) % Построение графика: подобранная кривая, f(x2) и экспериментальные точки. plot(x1, y1, ’-r’, x, y, ’*b’, x2, y2, ’sk’); % Вычисление индекса корреляции. R=sqrt(1-sum((y-yr ).^ 2)/sum((y-mean(y)).^2)) % Результаты вычислений c = 0.33503 0.90183 -0.69337 sum1=0.090533 x2= 2 3 4 y2= 0.65272 0.47033 0.30475 R= 0.99533
Таким образом подобрана зависимость $$Y = 0.33503x^{0.90183}e^{-0.9337x}$$. Вычислено ожидаемое значение в точках $$2, 3, 4: Y (2) = 0.65272, Y (3) = 0.47033, Y (4) = 0.30475$$. График подобранной зависимости вместе с экспериментальными точками и расчётными значениями изображён на рис. 11.3. Индекс корреляции равен 0.99533.
Способ 2.
Теперь рассмотрим решение задачи 11.2 путём решения системы 11.16. Решение с комментариями приведено в листинге 11.3. Результаты и графики при решении обоими способами полностью совпадают.
function s=f_mnk(c) % Переменные x,y являются глобальными, используются в нескольких функциях global x; global y; s =0; for i =1: length(x) s=s+(log(y(i))-c(1)-c(2)*log(x(i))-c(3)*x(i))^2; end end %___________________________________ global x; global y; % Определение координат экспериментальных точек x=[1 1.4 1.8 2.2 2.6 3 3.4 3.8 4.2 4.6 5 5.4 5.8]; y=[0.7 0.75 0.67 0.62 0.51 0.45 0.4 0.32 0.28 0.25 0.22 0.16 0.1]; % Формирование СЛАУ (11.16) G=[length(x) sum(log(x)) sum(x); sum(log(x)) sum(log(x).*log(x ))sum(x.*log(x)); sum(x) sum(x.*log(x)) sum(x.*x)]; H=[sum(log(y)); sum(log(y).*log(x)); sum(log(y).*x)]; % Решение СЛАУ методом Гаусса с помощью функции rref. C=rref([G H]); n=size(C); c=C(:, n(2)) % Вычисление суммарной квадратичной ошибки для подобранной % зависимости и вывод её на экран. sum1=f_mnk(c) % Формирование точек для построения графика подобранной кривой. x1 = 1:0.1:6; y1= exp(c(1)).*x1.^c(2).*exp(c(3).*x1); % Вычисление значений на подобранной кривой в заданных точках. yr=exp(c(1)).*x.^c(2).*exp(c(3).*x); % Вычисление ожидаемого значения подобранной функции в точках x=[2,3,4] x2 =[2 3 4 ] y2= exp(c(1)).*x2.^c(2).*exp(c(3).*x2) % Построение графика: подобранная кривая, f(x2) и экспериментальные точки. plot(x1, y1, ’-r’, x, y, ’*b’, x2, y2, ’sk’); % Вычисление индекса корреляции. R=sqrt(1-sum((y-yr).^2)/sum((y-mean(y)).^2))
Пример 11.3. В результате эксперимента получена табличная зависимость $$y(x)$$ (см. табл. 11.5). Подобрать аналитические зависимости $$f(x)=b_1+b_2x+b_3x^2+b_4x^3+b_5x^4+b_6x^5,g(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^5$$ и $$\varphi(x)=c_1+c_2x+c_4x^3+c_5x^5$$ методом наименьших квадратов. Пользуясь значением индекса корреляции выбрать наилучшую из них, с помощью которой вычислить ожидаемое значение в точках 1, 2.5, 4.8. Построить графики экспериментальных точек, подобранных зависимостей. На графиках отобразить рассчитанные значения в точках 1, 2.5, 4.8.
Как рассматривалось ранее, решать задачу подбора параметров полинома методом наименьших квадратов в Octave можно тремя способами.
(рис 11.3) График к примеру 11.2: экспериментальные точки и подобранная методом наименьших квадратов зависимость
| x | -2 | -1,3 | -0,6 | 0,1 | 0,8 | 1,5 | 2,2 | 2,9 | 3,6 | 4,3 | 5 | 5,7 | 6,4 |
| y | -10 | -5 | 0 | 0,7 | 0,8 | 2 | 3 | 5 | 8 | 30 | 60 | 100 | 238 |
Чтобы продемонстрировать использование всех трёх методов для подбора $$f(x)=\sum _{i=1}^6b_ix^{i-1}$$ воспользуемся функцией $$polyfit$$, для формирования коэффициентов функции $$g(x)$$ сформируем и решим систему уравнений (11.3), а функцию $$\varphi(x)$$ будем искать с помощью функции $$sqp$$.
Для формирования подбора коэффициентов функции $$g(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^5$$ сформируем систему уравнений. Составим функцию $$S(a_1,a_2,a_3,a_4,a_5)=\sum_{i=1}^n(y_i-a_1-a_2x_i-a_3x_i^2-a_4x_i^3-a_5x_i^5)^2$$. После дифференцирования $$S$$ по $$a_1,a_2,a_3,a_4$$ и $$a_5$$ система линейных алгебраических уравнений для вычисления параметров $$a_1,a_2,a_3,a_4,a_5$$ примет вид:
$$\left\{ \begin{aligned} a_1n+a_2\sum_{i=1}^nx_i+a_3\sum _{i=1}^nx_i^2+a_4\sum_{i=1}^nx_i^3+a_5\sum_{i=1}^nx_i^5=\sum_{i=1}^ny_i\\ a_1\sum_{i=1}^nx_i+a_2\sum_{i=1}^nx_i^2+a_3\sum_{i=1}^nx_i^3+a_4\sum_{i=1}^nx_i^4+a_5\sum_{i=1}^nx_i^6=\sum_{i=1}^ny_ix_i\\ a_1\sum_{i=1}^nx_i^2+a_2\sum_{i=1}^nx_i^3+a_3\sum_{i=1}^nx_i^4+a_4\sum_{i=1}^nx_i^5+a_5\sum_{i=1}^nx_i^7=\sum_{i=1}^ny_ix_i^2\\ a_1\sum_{i=1}^nx_i^3+a_2\sum_{i=1}^nx_i^4+a_3\sum_{i=1}^nx_i^5+a_4\sum_{i=1}^nx_i^6+a_5\sum_{i=1}^nx_i^8=\sum_{i=1}^ny_ix_i^3\\ a_1\sum_{i=1}^nx_i^5+a_2\sum_{i=1}^nx_i^6+a_3\sum_{i=1}^nx_i^7+a_4\sum_{i=1}^nx_i^8+a_5\sum_{i=1}^nx_i^{10}=\sum_{i=1}^ny_ix_i^5 \end{aligned} \right.$$Решив систему (11.21), найдём коэффициенты$$a_1,a_2,a_3,a_4$$ и $$a_5$$ функции$$g(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^5$$.
Для поиска функциональной зависимости вида $$\varphi(x)=c_1+c_2x+c_4x^3+c_5x^5$$ необходимо будет найти такие значения $$c_1,c_2,c_3,c_4$$ которых функция
$$S(c_1,c_2,c_3,c_4)=\sum_{i=1}^n\left(y_i-c_1-c_2x_i-c_3x_i^3-c_4x_i^5\right)^2$$принимала бы наименьшее значение.
После вывода необходимых формул приступим к реализации в Octave. Текст программы с очень подробными комментариями приведён в листинге 11.4.
% Функция для подбора зависимости fi(x) методом наименьших квадратов. function s=f_mnk(c) % Переменные x, y являются глобальными, используются в % функции f_mnk и главной функции. global x; global y; % Формирование суммы квадратов отклонений (11.22). s =0; for i =1: length(x) s=s+(y(i)-c(1)-c(2)*x(i)-c(3)*x(i)^3_c(4)*x(i)^5)^2; end end % Главная функция ———————————— global x; global y; % Определение координат экспериментальных точек x=[-2 -1.3 -0.6 0.1 0.8 1.5 2.2 2.9 3.6 4.3 5 5.7 6.4]; y=[-10 -5 0 0.7 0.8 2 3 5 8 30 60 100 2 3 8]; z =[1 2.5 4.8] % Подбор коэффициентов зависимости f(x) (полинома пятой степени) % методом наименьших квадратов, используя функцию polyfit. % Коэффициенты полинома будут хранится в переменной B. B=polyfit(x, y, 5) % Формирование точек для построения графиков подобранных функций. X1= -2:0.1:6.5; % Вычисление ординат точек графика первой функции f(x). Y1=polyval (B, X1); % Формирование системы (11.21) для подбора функции g(x). Здесь GGL — % матрица коэффициентов, H — вектор правых частей системы (11.21), % G — первые 4 строки и 4 столбца матрицы коэффициентов, G1 — пятый % столбец матрицы коэффициентов, G2 — пятая строка матрицы коэфф–тов. for i = 1:4 for j =1:4 G(i, j)= sum(x.^(i+j-2)); endfor endfor for i = 1:4 G1(i)= sum(x.^(i+5)); H(i)= sum(y.*x.^(i-1)); endfor for i =1:4 G2(i)= sum(x.^(i+4)); endfor G2(5)= sum(x.^10); % Формирование матрицы коэффициентов системы (11.21) из матриц % G, G1 и G2. GGL=[G G1’; G2] H(5)= sum(y.*x.^5); % Решение системы (11.21) методом обратной матрицы и % формирование коэффициентов А функции g(x). A=inv(GGL)*H’ % Подбор коэффициентов зависимости fi(x) методом наименьших квадратов, % используя функцию sqp. Коэффициенты функции будут хранится в перемен-% ной C. Задание начального значения вектора С, при неправильном его опре-% делении, экстремум функции может быть найден неправильно. C = [2; 1; 3; 1]; % Поиск вектора С, при котором функция (11.22) достигает своего % минимального значения, вектор С — коэффициенты функции fi. C=sqp (C, @f_mnk) % Вычисление ординат точек графика второй функции g(x). Y2=A(1)+A(2)*X1+A(3)*X1.^2+A(4)*X1.^3+A(5)*X1.^5; % Вычисление ординат точек графика третьей функции fi(x). Y3=C(1)+C(2)*X1+C(3)*X1.^3+C(4)*X1.^5; % Вычисление значений первой функции f(x) в заданных точках. yr1=polyval(B, x); % Вычисление значений второй функции g(x) в заданных точках. yr2=A(1)+A(2)*x+A(3)*x.^2+A(4)*x.^3+A(5)*x.^5; % Вычисление значений третьей функции fi(x) в заданных точках. yr3=C(1)+C(2)*x+C(3)*x.^3+C(4)*x.^5; % Вычисление индекса корреляции для первой функции f(x). R1=sqrt(1-sum((y-yr1).^2)/sum((y-mean(y)).^2)) % Вычисление индекса корреляции для второй функции g(x). R2=sqrt(1-sum((y-yr2).^2)/sum((y-mean(y)).^2)) % Вычисление индекса корреляции для третьей функции fi(x). R3=sqrt(1-sum((y-yr3).^2)/sum((y-mean(y)).^2)) % Сравнивая значения трёх индексов корреляции, выбираем наилучшую % функцию и с её помощью вычисляем ожидаемое значение в точках 1, 2.5, 4.8. if R1>R2 R1>R3 yz=polyval(B, z) "R1="; R1 endif if R2>R1 R2>R3 yz=C2(1)+C2(2)*z+C2(3)*z.^2+C2(4)*z.^3+C2(5)*z.^5 "R2="; R2 endif if R3>R1 R3>R2 yz=C(1)+C(2)*z+C(3)*z.^3+C(4)*z.^5 "R3="; R3 endif % Построение графика. plot(x, y, "*r;experiment;",X1, Y1, ’-b;f(x);’,X1, Y2, ’dr;g(x);’,X1 , Y3, ’ok;fi(x);’, z, yz, ’sb;f(z);’); grid(); % Результаты работы программы. z = 1.0000 2.5000 4.8000 B = 0.083039 -0.567892 0.906779 1.609432 -1.115925 -1.355075 GGL = 1.3000e+01 2.8600e+01 1.5210e+02 7.2701e+02 1.2793e+05 2.8600e+01 1.5210e+02 7.2701e+02 3.9868e+03 7.5030e+05 1.5210e+02 7.2701e+02 3.9868e+03 2.2183e+04 4.4706e+06 7.2701e+02 3.9868e+03 2.2183e+04 1.2793e+05 2.6938e+07 2.2183e+04 1.2793e+05 7.5030e+05 4.4706e+06 1.6383e+08 A = 9.4262e+00 -3.6516e+00 -5.7767e+00 1.7888e+00 -5.8179e-05 C = -1.030345 5.080391 -0.609721 0.033534 R1 = 0.99690 R2 = 0.98136 R3 = 0.99573 yz = -0.43964 6.008544 0.77972
На рисунке 11.4 представлено графическое решение задачи.
(рис 11.4) Графическое решение к примеру 11.3
Рассмотренная задача демонстрирует основные приёмы подбора зависимости методом наименьших квадратов. Авторы рекомендует внимательно рассмотреть её для понимания методов решения подобных задач в Octave.
В заключении авторы позволят несколько советов по решению задачи аппроксимации.
Метод наименьших квадратов (МНК) позволяет по экспериментальным данным подобрать такую аналитическую функцию, которая проходит настолько близко к экспериментальным точкам, насколько это возможно.
В общем случае задачу можно сформулировать следующим образом. Пусть в результате эксперимента были получены некая экспериментальная зависимость $$y(x)$$, представленная в таблице 11.1. Необходимо построить аналитическую зависимость $$f(x,a_{1},a_{2},\dots, a_{k})$$, наиболее точно описывающую результаты эксперимента. Для построения параметров функции $$f(x,a_{1},a_{2},\dots, a_{k})$$ будем использовать метод наименьших квадратов. Идея метода наименьших квадратов заключается в том, что функцию $$f(x,a_{1},a_{2},\dots, a_{k})$$ необходимо подобрать таким образом, чтобы сумма квадратов отклонений измеренных значений $$Y_i=f(x,a_{1},a_{2},\dots, a_{k})$$ была бы наименьшей (см. рис. 11.1):
| x | $$x_1$$ | $$x_2$$ | $$x_3$$ | ... | $$x_n-1$$ | $$x_n$$ |
| y | $$y_1$$ | $$y_2$$ | $$y_3$$ | ... | $$y_n-1$$ | $$y_n$$ |
Задача состоит из двух этапов:
Математически задача подбора коэффициентов зависимости сводится к определению коэффициентов $$a_i$$ из условия (11.1). В Octave её можно решать несколькими способами:
Вспомним некоторые сведения из высшей математики, необходимые для решения задачи подбора зависимости методом наименьших квадратов.
Достаточным условием минимума функции $$S(a_{1},a_{2},\dots, a_{k})$$ является равенство нулю всех её частных производных. Поэтому задача поиска минимума функции (11.1) эквивалентна решению системы
(рис 11.1) Геометрическая интерпретация МНК
алгебраических уравнений:
$$\left\{ \begin{matrix} \frac{\partial S}{\partial a_1}=0\\ \frac{\partial S}{\partial a_2}=0\\ \dots\\ \frac{\partial S}{\partial a_k}=0 \end{matrix} \right.$$Если параметры $$a_i$$ входят в зависимость $$Y=f(x,a_{1},a_{2},\dots, a_{k})$$ линейно, то получим систему (11.3) из k линейных уравнений с $$k$$ неизвестными.
$$\left\{ \begin{matrix} \displaystyle \sum_{i=1}^n2\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)\frac{\partial f}{\partial a_1}=0\\ \displaystyle \sum_{i=1}^n2\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)\frac{\partial f}{\partial a_2}=0\\ \hdotsfor{1}\\ \displaystyle \sum_{i=1}^n2\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)\frac{\partial f}{\partial a_k}=0 \end{matrix} \right.$$Составим систему (11.3) для наиболее часто используемых функций.
Для подбора параметров линейной функции $$Y=a_1+a_2x$$, составим функцию (11.1) для линейной зависимости:
$$S(a_1,a_2)=\sum_{i=1}^n\left(y_i-a_1-a_2x_i\right)^2\to \min.$$Продифференцировав функцию $$S$$ по $$a_1$$ и $$a_2$$, получим систему уравнений:
$$\left\{ \begin{matrix} \displaystyle 2\sum_{i=1}^n\left(y_i-a_1-a_2x_i\right)(-1)=0\\ \displaystyle 2\sum_{i=1}^n\left(y_i-a_1-a_2x_i\right)(-x_i)=0 \end{matrix} \right. \Rightarrow \left\{ \begin{matrix} \displaystyle a_1n+a_2\sum_{i=1}^nx_i=\sum_{i=1}^ny_i\\ \displaystyle a_1\sum_{i=1}^nx_i+a_2\sum_{i=1}^nx_i^2=\sum_{i=1}^ny_ix_i \end{matrix} \right.$$решив которую, определим коэффициенты функции $$Y=a_1+a_2x$$:
$$\left\{ \begin{matrix} \displaystyle a_1=\displaystyle \frac{\displaystyle\sum_{i=1}^ny_i}n-a_2\displaystyle \frac{\displaystyle \sum_{i=1}^nx_i}n\\ \displaystyle a_2=\displaystyle \frac{n\displaystyle\sum_{i=1}^ny_ix_i-\sum_{i=1}^ny_i\sum_{i=1}^nx_i}{n\displaystyle\sum_{i=1}^nx_i^2-\left(\sum_{i=1}^nx_i\right)^2} \end{matrix} \right.$$Для определения параметров зависимости $$Y=a_1+a_2x+a_3x^2$$ составим функцию $$S(a_1,a_2,a_3)$$(11.1):
$$S(a_1,a_2,a_3)=\sum_{i=1}^n\left(y_i-a_1-a_2x_i-a_3x_i^2\right)^2\to \min$$После дифференцирования S по $$a_i$$,$$a_2$$ и $$a_3$$ получим систему линейных алгебраических уравнений:
$$\left\{ \begin{matrix} \displaystyle a_{1}n+a_{2}\sum_{i=1}^{n}x_{i}+a_{3}\sum_{i=1}^{n}x_{i}^{2}=\sum_{i=1}^{n}y_{i}\\ \displaystyle a_{1}\sum_{i=1}^{n}x_{i}+a_{2}\sum_{i=1}^{n}x_{i}^{2}+a_{3}\sum_{i=1}^{n}x_{i}^{3}=\sum_{i=1}^{n}y_{i}x_{i}\\ \displaystyle a_{1}\sum_{i=1}^{n}x_{i}^{2}+a_{2}\sum_{i=1}^{n}x_{i}^{3}+a_{3}\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}y_{i}x_{i}^{2} \end{matrix} \right.$$Решив систему (11.8), найдём значения параметров a$$a_i$$,$$a_2$$ и $$a_3$$.
Аналогично определим параметры многочлена третьей степени:$$S(a_1,a_2,a_3,a_4)$$. Составим функцию $$S(a_1,a_2,a_3,a_4)$$:
$$S(a_1,a_2,a_3,a_4)=\sum_{i=1}^n\left(y_i-a_1-a_2x_i-a_3x_i^2-a_4x_i^3\right)^2\to\min$$После дифференцирования $$S$$ по $$a_i$$,$$a_2$$,$$a_3$$ и $$a_4$$, система линейных алгебраических уравнений для вычисления параметров a $$a_1,a_2,a_3,a_4$$ примет вид:
$$\left\{ \begin{matrix} \displaystyle a_{1}n+a_{2}\sum_{i=1}^{n}x_{i}+a_{3}\sum_{i=1}^{n}x_{i}^{2}+a_{4}\sum_{i=1}^{n}x_{i}^{3}=\sum_{i=1}^{n}y_{i}\\ \displaystyle a_{1}\sum _{i=1}^{n}x_{i}+a_{2}\sum _{i=1}^{n}x_{i}^{2}+a_{3}\sum_{i=1}^{n}x_{i}^{3}+a_{4}\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}y_{i}x_{i}\\ \displaystyle a_{1}\sum _{i=1}^{n}x_{i}^{2}+a_{2}\sum_{i=1}^{n}x_{i}^{3}+a_{3}\sum _{i=1}^{n}x_{i}^{4}+a_{4}\sum_{i=1}^{n}x_{i}^{5}=\sum_{i=1}^{n}y_{i}x_{i}^{2}\\ \displaystyle a_{1}\sum_{i=1}^{n}x_{i}^{3}+a_{2}\sum_{i=1}^{n}x_{i}^{4}+a_{3}\sum_{i=1}^{n}x_{i}^{5}+a_{4}\sum_{i=1}^{n}x_{i}^{6}=\sum_{i=1}^{n}y_{i}x_{i}^{3} \end{matrix} \right.$$Решив систему (11.10), найдём коэффициенты $$a_i$$,$$a_2$$,$$a_3$$ и $$a_4$$.
В общем случае система уравнений для вычисления параметров $$a_i$$ многочлена k-й степени $$Y=\sum_{i=1}^{k+1}a_{i}x^{i-1}$$ имеет вид:
$$\left\{ \begin{matrix} \displaystyle a_1n+a_2\sum_{i=1}^nx_i+a_3\sum_{i=1}^nx_i^2+\ldots +a_{k+1}\sum_{i=1}^nx_i^k=\sum_{i=1}^ny_i\\ \displaystyle a_1\sum_{i=1}^nx_i+a_2\sum_{i=1}^nx_i^2+a_3\sum_{i=1}^nx_i^3+\ldots +a_{k+1}\sum_{i=1}^nx_i^{k+1}=\sum_{i=1}^ny_ix_i\\ \hdotsfor{1}\\ \displaystyle a_1\sum_{i=1}^nx_i^{k+1}+a_2\sum_{i=1}^nx_i^{k+2}+a_3\sum_{i=1}^nx_i^{k+3}+\ldots +a_{k+1}\sum_{i=1}^nx_i^{2k}=\sum_{i=1}^ny_ix_i^k \end{matrix} \right.$$В матричном виде систему (11.11) можно записать
$$Ca = g,$$Элементы матрицы C и вектора g рассчитываются по формулам
$$C_{i,j}=\sum_{k=1}^nx_k^{i+j-2},\quad i=1,\dots,k+1,\ j=1,\dots,k+1$$ $$g_i=\sum_{k=1}^ny_kx_k^{i-1},\quad i=1,\dots,k+1$$Решив систему (11.12), определим параметры зависимости $$Y=a_1+a_2x+a_3x^2+\dots+a_{k+1}x^k$$.
$$Y=ax^be^{cx}$$
Параметры $$b$$ и $$c$$ входят в зависимость $$Y=ax^be^{cx}$$ нелинейным образом. Чтобы избавиться от нелинейности предварительно
Составим функцию $$S(A, b, c)$$ по формуле (11.1):
$$S(A,b,c)=\sum_{i=1}^n\left(Y1_i-A-b\ln x_i-cx_i\right)^2\to\min$$После дифференцирования получим систему трёх линейных алгебраических уравнений для определения коэффициентов $$A, b, c$$.
$$\left\{ \begin{matrix} \displaystyle nA+b\sum_{i=1}^n\ln x_i+c\sum_{i=1}^nx_i=\sum_{i=1}^nY1_i\\ \displaystyle A\sum_{i=1}^n\ln x_i+b\sum_{i=1}^n(\ln x_i)^2+b\sum_{i=1}^nx_i\ln x_i=\sum_{i=1}^nY1_i\ln x_i\\ \displaystyle A\sum_{i=1}^nx_i+b\sum_{i=1}^nx_i\ln x_i+c\sum_{i=1}^nx_i^2=\sum_{i=1}^nY1_ix_i \end{matrix} \right.$$После решения системы (11.16) необходимо вычислить значение коэффициента a по формуле $$a=e^A$$.
Для вычисления параметров функции $$Y = ax^b$$ необходимо предварительно её прологарифмировать $$\ln Y=\ln ax^b=\ln a+b\ln x$$. После чего замена $$Z = lnY, X = lnx, A = lna$$ приводит заданную функцию к линейному виду $$Z = bX + A$$, где коэффициенты $$A$$ и $$b$$ вычисляются по формулам (11.6) и, соответственно, $$a=e^A$$.
Аналогично можно подобрать параметры функции вида $$Y = ae^{bx}$$. Прологарифмируем заданную функцию $$\ln y=\ln a+bx\ln e,\ln y=\ln a+bx$$. Проведём замену $$Y = lny, A = lna$$ и получим линейную зависимость $$Y = bx + A$$. По формулам (11.6) найдём $$A$$ и $$b$$, а затем вычислим $$a=e^A$$.
Рассмотрим ещё ряд зависимостей, которые сводятся к линейной.
Для подбора параметров функции $$Y=\frac{1}{ax+b}$$ сделаем замену $$Y=\frac{1}{Y}$$. В результате получим линейную зависимость Z = ax + b. Функция $$Y=\frac{x}{ax+b}$$ заменами $$Z=\frac{1}{Y}$$, $$X=\frac{1}{x}$$ сводится к линейной Z = a + bX. Для определения коэффициентов функциональной зависимости $$Y=\frac{1}{ae^{-x}+b}$$ необходимо сделать следующие замены $$Z=\frac{1}{Y},X=e^{-x}$$. В результате также получим линейную функцию Z = aX + b.
Аналогичными приёмами (логарифмированием, заменами и т. п.) можно многие подбираемые зависимости преобразовать к такому виду, что получаемая при решении задачи оптимизации система (11.2) была системой линейных алгебраических уравнений. При использовании Octave можно напрямую решать задачу подбора параметров, как задачу оптимизации (11.1) с использованием функции sqp.
После нахождения параметров зависимости $$f(x,a_{1},a_{2},\dots, a_{k})$$ возникает вопрос насколько адекватно описывает подобранная зависимость экспериментальные данные. Чем ближе величина
$$S=\sum_{i=1}^n\left(y_i-f(x_i,a_1,a_2,\dots,a_k)\right)^2$$называемая суммарной квадратичной ошибкой, к нулю, тем точнее подобранная кривая описывает экспериментальные данные.
Линия, описываемая уравнением вида $$y=a_{1}+a_{2}x$$, называется линией регрессии $$y$$ на $$x$$, параметры $$a_1$$ и $$a_2$$ называются коэффициентами регрессии и определяются формулами (11.6).
Чем меньше величина $$S=\sum _{i=1}^{n}(y_{i}-a_{1}-a_{2}x_{i})^{2}$$, тем более обоснованно предположение, что экспериментальные данные описываются линейной функцией. Существует показатель, характеризующий тесноту линейной связи между $$x$$ и $$y$$, который называется коэффициентом корреляции и рассчитывается по формуле:
$$r=\frac{\displaystyle \sum_{i=1}^n\left(x_i-M_x\right)\left(y_i-M_y\right)} {\sqrt{\displaystyle\sum_{i=1}^n\left(x_i-M_x\right)^2\sum_{i=1}^n\left(y_i-M_y\right)^2}}, \ M_{x}=\frac{\displaystyle\sum_{i=1}^nx_i}{n},\ M_{y}=\frac{\displaystyle\sum_{i=1}^ny_i}{n}$$Значение коэффициента корреляции удовлетворяет соотношению $$-1\le r\le 1$$.
Чем меньше отличается абсолютная величина $$r$$ от единицы, тем ближе к линии регрессии располагаются экспериментальные точки. Если $$|r| = 1$$, то все экспериментальные точки находятся на линии регрессии. Если коэффициент корреляции близок к нулю, то это означает, что между $$x$$ и $$y$$ не существует линейной связи, но между ними может существовать зависимость, отличная от линейной.
Для того, чтобы проверить, значимо ли отличается от нуля коэффициент корреляции, можно использовать критерий Стьюдента. Вычисленное значение критерия определяется по формуле:
$$t=r\sqrt{\frac{n-2}{1-r^2}}$$Рассчитанное по формуле (11.19) значение $$t$$ сравнивается со значением, взятым из таблицы распределения Стьюдента (см. табл. 11.2) в соответствии с уровнем значимости $$p$$ (стандартное значение $$p = 0.95$$) и числом степеней свободы $$k = n - 2$$. Если полученная по формуле (11.19) величина $$t$$ больше табличного значения, то коэффициент корреляции значимо отличен от нуля.
| k\p | 0,99 | 0,98 | 0,95 | 0,90 | 0,80 | 0,70 | 0,60 |
| 1 | 63,657 | 31,821 | 12,706 | 6,314 | 3,078 | 1,963 | 1,376 |
| 2 | 9,925 | 6,965 | 4,303 | 2,920 | 1,886 | 1,386 | 1,061 |
| 3 | 5,841 | 4,541 | 3,182 | 2,353 | 1,638 | 1,250 | 0,978 |
| 4 | 4,604 | 3,747 | 2,776 | 2,132 | 1,533 | 1,190 | 0,941 |
| 5 | 4,032 | 3,365 | 2,571 | 2,05 | 1,476 | 1,156 | 0,920 |
| 6 | 3,707 | 3,141 | 2,447 | 1,943 | 1,440 | 1,134 | 0,906 |
| 7 | 3,499 | 2,998 | 2,365 | 1,895 | 1,415 | 1,119 | 0,896 |
| 8 | 3,355 | 2,896 | 2,896 | 1,860 | 1,387 | 1,108 | 0,889 |
| 9 | 3,250 | 2,821 | 2,261 | 1,833 | 1,383 | 1,100 | 0,883 |
| 10 | 3,169 | 2,764 | 2,228 | 1,812 | 1,372 | 1,093 | 0,879 |
| 11 | 3,106 | 2,718 | 2,201 | 1,796 | 1,363 | 1,088 | 0,876 |
| 12 | 3,055 | 2,681 | 2,179 | 1,782 | 1,356 | 1,083 | 0,873 |
| 13 | 3,012 | 2,681 | 2,179 | 1,782 | 1,356 | 1,083 | 0,873 |
| 14 | 2,977 | 2,624 | 2,145 | 1,761 | 1,345 | 1,076 | 0,868 |
| 15 | 2,947 | 2,602 | 2,131 | 1,753 | 1,341 | 1,074 | 0,866 |
| 16 | 2,921 | 2,583 | 2,120 | 1,746 | 1,337 | 1,071 | 0,865 |
| 17 | 2,898 | 2,567 | 2,110 | 1,740 | 1,333 | 1,069 | 0,863 |
| 18 | 2,878 | 2,552 | 2,101 | 1,734 | 1,330 | 1,067 | 0,862 |
| 19 | 2,861 | 2,539 | 2,093 | 1,729 | 1,328 | 1,066 | 0,861 |
| 20 | 2,845 | 2,528 | 2,086 | 1,725 | 1,325 | 1,064 | 0,860 |
| 21 | 2,831 | 2,518 | 2,080 | 1,721 | 1,323 | 1,063 | 0,859 |
| 22 | 2,819 | 2,508 | 2,074 | 1,717 | 1,321 | 1,061 | 0,858 |
| 23 | 2,807 | 2,500 | 2,069 | 1,714 | 1,319 | 1,060 | 0,858 |
| 24 | 2,797 | 2,492 | 2,064 | 1,711 | 1,318 | 1,059 | 0,857 |
| 25 | 2,779 | 2,485 | 2,060 | 1,708 | 1,316 | 1,058 | 0,856 |
| 26 | 2,771 | 2,479 | 2,056 | 1,706 | 1,315 | 1,058 | 0,856 |
| 27 | 2,763 | 2,473 | 2,052 | 1,703 | 1,314 | 1,057 | 0,855 |
| 28 | 2,756 | 2,467 | 2,048 | 1,701 | 1,313 | 1,056 | 0,855 |
| 29 | 2,750 | 2,462 | 2,045 | 1,699 | 1,311 | 1,055 | 0,854 |
| 30 | 2,704 | 2,457 | 2,042 | 1,697 | 1,310 | 1,055 | 0,854 |
| 40 | 2,660 | 2,423 | 2,021 | 1,684 | 1,303 | 1,050 | 0,851 |
| 60 | 2,612 | 2,390 | 2,000 | 1,671 | 1,296 | 1,046 | 0,848 |
| 120 | 2,617 | 2,358 | 1,980 | 1,980 | 1,289 | 1,041 | 0,845 |
| $$\infty$$ | 2,576 | 2,326 | 1,960 | 1,645 | 1,282 | 1,036 | 0,842 |
Коэффициент корреляции $$r$$ применяется только в тех случаях, когда между данными существует прямолинейная связь. Если же связь нелинейная, то для выявления тесноты связи между переменными $$y$$ и $$x$$ в случае нелинейной зависимости пользуются индекс корреляции. Он показывает тесноту связи между фактором $$x$$ и зависимой переменной $$y$$ и рассчитывается по формуле:
$$R=\sqrt{1-\frac{\displaystyle \sum _{i=1}^{n}\left(y_{i}-Y_{i}\right)^{2}}{\sum_{i=1}^{n}\left(y_{i}-M_{y}\right)^{2}}}$$где $$y$$ — экспериментальные значения, $$Y $$— теоретические значения (рассчитанные по подобранной методом наименьших квадратов формуле), $$M_y$$ — среднее значение $$y$$.
Индекс корреляции лежит в пределах от 0 до 1. При наличии функциональной зависимости индекс корреляции близок к 1. При отсутствии связи $$R$$ практически равен нулю. Если коэффициент корреляции r является мерой тесноты связи только для линейной формы связи, то индекс корреляции $$R$$ — как для линейной, так и для нелинейной. При прямолинейной связи коэффициент корреляции по своей абсолютной величине равен индексу корреляции: $$|r| = R$$.
| t | $$0^o$$ | $$4^o$$ | $$10^o$$ | $$15^o$$ | $$21^o$$ | $$29^o$$ | $$36^o$$ | $$51^o$$ | $$68^o$$ |
| P | -66,7 | -71,0 | 76,3 | 80,6 | 85,7 | 92,9 | 99,4 | 113,6 | 125,1 |
Для решения задач подбора аналитических зависимостей по экспериментальным данным можно использовать следующие функции Octave:
Пример 11.1. В ->Основах химии-> Д.И. Менделеева приводятся данные о растворимости азотнокислого натрия $$NaNO_3$$ в зависимости от температуры воды. В 100 частях воды (табл. 11.3) растворяется следующее число условных частей $$NaNO_3$$ при соответствующих температурах. Требуется определить растворимость азотнокислого натрия при температуре $$t = 32°C$$ в случае линейной зависимости и найти коэффициент корреляции.
Решение задачи 11.1 с комментариями приведено в листинге 11.1.
(рис 11.2) Иллюстрация к примеру 11.1
% Ввод экспериментальных данных X=[0 4 10 15 21 29 36 51 68]; Y=[66.7 71.0 76.3 80.6 85.7 92.9 99. 4 113.6 125.1]; % Вычисление вектора коэффициентов полинома y=a1*x+a2 [ a]= polyfit (X, Y, 1) % Вычисление значения полинома y=a1*x+a2 в точке t=32 t =32; yt=a(1)*t+a(2) % Построение графика полинома y=a1*x+a2, экспериментальных точек % и значения в заданной точке в одной графической области x =0:68; y=a(1)*x+a(2); plot (X, Y, ’ok’, x, y, ’-k’, t, yt, ’*k’) grid on; % Вычисление коэффициента корреляции k =cor(x, y) % Результаты вычислений a = 0.87064 67.50779 yt = 95.368 k = 1
На рис. 11.2 приведено графическое решение этой задачи, изображены экспериментальные точки, линия регрессии $$y=a_1x+a_2$$, на котором отмечена точка $$t = 32$$.
Пример 11.2. В результате эксперимента получена табличная зависимость $$y(x)$$ (см. табл. 11.4). Подобрать аналитическую зависимость $$Y=ax^be^{cx}$$ методом наименьших квадратов. Вычислить ожидаемое значение в точках 2, 3, 4. Вычислить индекс корреляции. Решение задачи подбора параметров функции $$f(x)=ax^be^{cx}$$ в Octave возможно двумя способами:
| x | 1 | 1,4 | 1,8 | 2,2 | 2,6 | 3 | 3,4 | 3,8 | 4,2 | 4,6 | 5 | 5,4 | 5,8 |
| y | 0,7 | 0,75 | 0,67 | 0,62 | 0,51 | 0,45 | 0,4 | 0,32 | 0,28 | 0,25 | 0,22 | 0,16 | 0,1 |
Рассмотрим последовательно оба варианта решения задачи.
Способ 1.
Функция (11.15) реализована в Octave с помощью функции $$f\_mnk$$. Полный текст программы решения задачи способом 1 с комментариями приведён на листинге 11.2. Вместо коэффициентов $$A, b, c$$ из формул (11.15)–(11.16) в программе на Octave используется массив $$c$$.
function s=f_mnk( c ) % Переменные x,y являются глобальными, используются в нескольких функциях global x; global y; s =0; for i =1: length(x) s=s +(log(y(i))-c(1)-c(2)*log(x(i))*c(3)*x(i))^2; end end % _________________________________________________ global x; global y; % Задание начального значения вектора c, при неправильном его % определении, экстремум может быть найден неправильно. c = [2; 1; 3]; % Определение координат экспериментальных точек x=[1 1.4 1.8 2.2 2.6 3 3.4 3.8 4.2 4.6 5 5.4 5.8]; y=[0.7 0.75 0.67 0.62 0.51 0.45 0.4 0.32 0.28 0.25 0.22 0.16 0.1]; % Решение задачи оптимизации функции 11.15 с помощью sqp. c=sqp(c,@f_mnk) % Вычисление суммарной квадратичной ошибки для подобранной % зависимости и вывод её на экран. sum1=f_mnk( c ) % Формирование точек для построения графика подобранной кривой. x1 = 1:0.1:6; y1= exp(c(1)).*x1.^ c(2).*exp(c(3).*x1); % Вычисление значений на подобранной кривой в заданных точках. yr=exp(c(1)).*x.^ c(2).*exp(c(3).*x ); % Вычисление ожидаемого значения подобранной функции в точках x=[2,3,4] x2 =[2 3 4 ] y2= exp(c(1)).*x2.^c(2).*exp(c(3).*x2) % Построение графика: подобранная кривая, f(x2) и экспериментальные точки. plot(x1, y1, ’-r’, x, y, ’*b’, x2, y2, ’sk’); % Вычисление индекса корреляции. R=sqrt(1-sum((y-yr ).^ 2)/sum((y-mean(y)).^2)) % Результаты вычислений c = 0.33503 0.90183 -0.69337 sum1=0.090533 x2= 2 3 4 y2= 0.65272 0.47033 0.30475 R= 0.99533
Таким образом подобрана зависимость $$Y = 0.33503x^{0.90183}e^{-0.9337x}$$. Вычислено ожидаемое значение в точках $$2, 3, 4: Y (2) = 0.65272, Y (3) = 0.47033, Y (4) = 0.30475$$. График подобранной зависимости вместе с экспериментальными точками и расчётными значениями изображён на рис. 11.3. Индекс корреляции равен 0.99533.
Способ 2.
Теперь рассмотрим решение задачи 11.2 путём решения системы 11.16. Решение с комментариями приведено в листинге 11.3. Результаты и графики при решении обоими способами полностью совпадают.
function s=f_mnk(c) % Переменные x,y являются глобальными, используются в нескольких функциях global x; global y; s =0; for i =1: length(x) s=s+(log(y(i))-c(1)-c(2)*log(x(i))-c(3)*x(i))^2; end end %___________________________________ global x; global y; % Определение координат экспериментальных точек x=[1 1.4 1.8 2.2 2.6 3 3.4 3.8 4.2 4.6 5 5.4 5.8]; y=[0.7 0.75 0.67 0.62 0.51 0.45 0.4 0.32 0.28 0.25 0.22 0.16 0.1]; % Формирование СЛАУ (11.16) G=[length(x) sum(log(x)) sum(x); sum(log(x)) sum(log(x).*log(x ))sum(x.*log(x)); sum(x) sum(x.*log(x)) sum(x.*x)]; H=[sum(log(y)); sum(log(y).*log(x)); sum(log(y).*x)]; % Решение СЛАУ методом Гаусса с помощью функции rref. C=rref([G H]); n=size(C); c=C(:, n(2)) % Вычисление суммарной квадратичной ошибки для подобранной % зависимости и вывод её на экран. sum1=f_mnk(c) % Формирование точек для построения графика подобранной кривой. x1 = 1:0.1:6; y1= exp(c(1)).*x1.^c(2).*exp(c(3).*x1); % Вычисление значений на подобранной кривой в заданных точках. yr=exp(c(1)).*x.^c(2).*exp(c(3).*x); % Вычисление ожидаемого значения подобранной функции в точках x=[2,3,4] x2 =[2 3 4 ] y2= exp(c(1)).*x2.^c(2).*exp(c(3).*x2) % Построение графика: подобранная кривая, f(x2) и экспериментальные точки. plot(x1, y1, ’-r’, x, y, ’*b’, x2, y2, ’sk’); % Вычисление индекса корреляции. R=sqrt(1-sum((y-yr).^2)/sum((y-mean(y)).^2))
Пример 11.3. В результате эксперимента получена табличная зависимость $$y(x)$$ (см. табл. 11.5). Подобрать аналитические зависимости $$f(x)=b_1+b_2x+b_3x^2+b_4x^3+b_5x^4+b_6x^5,g(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^5$$ и $$\varphi(x)=c_1+c_2x+c_4x^3+c_5x^5$$ методом наименьших квадратов. Пользуясь значением индекса корреляции выбрать наилучшую из них, с помощью которой вычислить ожидаемое значение в точках 1, 2.5, 4.8. Построить графики экспериментальных точек, подобранных зависимостей. На графиках отобразить рассчитанные значения в точках 1, 2.5, 4.8.
Как рассматривалось ранее, решать задачу подбора параметров полинома методом наименьших квадратов в Octave можно тремя способами.
(рис 11.3) График к примеру 11.2: экспериментальные точки и подобранная методом наименьших квадратов зависимость
| x | -2 | -1,3 | -0,6 | 0,1 | 0,8 | 1,5 | 2,2 | 2,9 | 3,6 | 4,3 | 5 | 5,7 | 6,4 |
| y | -10 | -5 | 0 | 0,7 | 0,8 | 2 | 3 | 5 | 8 | 30 | 60 | 100 | 238 |
Чтобы продемонстрировать использование всех трёх методов для подбора $$f(x)=\sum _{i=1}^6b_ix^{i-1}$$ воспользуемся функцией $$polyfit$$, для формирования коэффициентов функции $$g(x)$$ сформируем и решим систему уравнений (11.3), а функцию $$\varphi(x)$$ будем искать с помощью функции $$sqp$$.
Для формирования подбора коэффициентов функции $$g(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^5$$ сформируем систему уравнений. Составим функцию $$S(a_1,a_2,a_3,a_4,a_5)=\sum_{i=1}^n(y_i-a_1-a_2x_i-a_3x_i^2-a_4x_i^3-a_5x_i^5)^2$$. После дифференцирования $$S$$ по $$a_1,a_2,a_3,a_4$$ и $$a_5$$ система линейных алгебраических уравнений для вычисления параметров $$a_1,a_2,a_3,a_4,a_5$$ примет вид:
$$\left\{ \begin{aligned} a_1n+a_2\sum_{i=1}^nx_i+a_3\sum _{i=1}^nx_i^2+a_4\sum_{i=1}^nx_i^3+a_5\sum_{i=1}^nx_i^5=\sum_{i=1}^ny_i\\ a_1\sum_{i=1}^nx_i+a_2\sum_{i=1}^nx_i^2+a_3\sum_{i=1}^nx_i^3+a_4\sum_{i=1}^nx_i^4+a_5\sum_{i=1}^nx_i^6=\sum_{i=1}^ny_ix_i\\ a_1\sum_{i=1}^nx_i^2+a_2\sum_{i=1}^nx_i^3+a_3\sum_{i=1}^nx_i^4+a_4\sum_{i=1}^nx_i^5+a_5\sum_{i=1}^nx_i^7=\sum_{i=1}^ny_ix_i^2\\ a_1\sum_{i=1}^nx_i^3+a_2\sum_{i=1}^nx_i^4+a_3\sum_{i=1}^nx_i^5+a_4\sum_{i=1}^nx_i^6+a_5\sum_{i=1}^nx_i^8=\sum_{i=1}^ny_ix_i^3\\ a_1\sum_{i=1}^nx_i^5+a_2\sum_{i=1}^nx_i^6+a_3\sum_{i=1}^nx_i^7+a_4\sum_{i=1}^nx_i^8+a_5\sum_{i=1}^nx_i^{10}=\sum_{i=1}^ny_ix_i^5 \end{aligned} \right.$$Решив систему (11.21), найдём коэффициенты$$a_1,a_2,a_3,a_4$$ и $$a_5$$ функции$$g(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^5$$.
Для поиска функциональной зависимости вида $$\varphi(x)=c_1+c_2x+c_4x^3+c_5x^5$$ необходимо будет найти такие значения $$c_1,c_2,c_3,c_4$$ которых функция
$$S(c_1,c_2,c_3,c_4)=\sum_{i=1}^n\left(y_i-c_1-c_2x_i-c_3x_i^3-c_4x_i^5\right)^2$$принимала бы наименьшее значение.
После вывода необходимых формул приступим к реализации в Octave. Текст программы с очень подробными комментариями приведён в листинге 11.4.
% Функция для подбора зависимости fi(x) методом наименьших квадратов. function s=f_mnk(c) % Переменные x, y являются глобальными, используются в % функции f_mnk и главной функции. global x; global y; % Формирование суммы квадратов отклонений (11.22). s =0; for i =1: length(x) s=s+(y(i)-c(1)-c(2)*x(i)-c(3)*x(i)^3_c(4)*x(i)^5)^2; end end % Главная функция ———————————— global x; global y; % Определение координат экспериментальных точек x=[-2 -1.3 -0.6 0.1 0.8 1.5 2.2 2.9 3.6 4.3 5 5.7 6.4]; y=[-10 -5 0 0.7 0.8 2 3 5 8 30 60 100 2 3 8]; z =[1 2.5 4.8] % Подбор коэффициентов зависимости f(x) (полинома пятой степени) % методом наименьших квадратов, используя функцию polyfit. % Коэффициенты полинома будут хранится в переменной B. B=polyfit(x, y, 5) % Формирование точек для построения графиков подобранных функций. X1= -2:0.1:6.5; % Вычисление ординат точек графика первой функции f(x). Y1=polyval (B, X1); % Формирование системы (11.21) для подбора функции g(x). Здесь GGL — % матрица коэффициентов, H — вектор правых частей системы (11.21), % G — первые 4 строки и 4 столбца матрицы коэффициентов, G1 — пятый % столбец матрицы коэффициентов, G2 — пятая строка матрицы коэфф–тов. for i = 1:4 for j =1:4 G(i, j)= sum(x.^(i+j-2)); endfor endfor for i = 1:4 G1(i)= sum(x.^(i+5)); H(i)= sum(y.*x.^(i-1)); endfor for i =1:4 G2(i)= sum(x.^(i+4)); endfor G2(5)= sum(x.^10); % Формирование матрицы коэффициентов системы (11.21) из матриц % G, G1 и G2. GGL=[G G1’; G2] H(5)= sum(y.*x.^5); % Решение системы (11.21) методом обратной матрицы и % формирование коэффициентов А функции g(x). A=inv(GGL)*H’ % Подбор коэффициентов зависимости fi(x) методом наименьших квадратов, % используя функцию sqp. Коэффициенты функции будут хранится в перемен-% ной C. Задание начального значения вектора С, при неправильном его опре-% делении, экстремум функции может быть найден неправильно. C = [2; 1; 3; 1]; % Поиск вектора С, при котором функция (11.22) достигает своего % минимального значения, вектор С — коэффициенты функции fi. C=sqp (C, @f_mnk) % Вычисление ординат точек графика второй функции g(x). Y2=A(1)+A(2)*X1+A(3)*X1.^2+A(4)*X1.^3+A(5)*X1.^5; % Вычисление ординат точек графика третьей функции fi(x). Y3=C(1)+C(2)*X1+C(3)*X1.^3+C(4)*X1.^5; % Вычисление значений первой функции f(x) в заданных точках. yr1=polyval(B, x); % Вычисление значений второй функции g(x) в заданных точках. yr2=A(1)+A(2)*x+A(3)*x.^2+A(4)*x.^3+A(5)*x.^5; % Вычисление значений третьей функции fi(x) в заданных точках. yr3=C(1)+C(2)*x+C(3)*x.^3+C(4)*x.^5; % Вычисление индекса корреляции для первой функции f(x). R1=sqrt(1-sum((y-yr1).^2)/sum((y-mean(y)).^2)) % Вычисление индекса корреляции для второй функции g(x). R2=sqrt(1-sum((y-yr2).^2)/sum((y-mean(y)).^2)) % Вычисление индекса корреляции для третьей функции fi(x). R3=sqrt(1-sum((y-yr3).^2)/sum((y-mean(y)).^2)) % Сравнивая значения трёх индексов корреляции, выбираем наилучшую % функцию и с её помощью вычисляем ожидаемое значение в точках 1, 2.5, 4.8. if R1>R2 R1>R3 yz=polyval(B, z) "R1="; R1 endif if R2>R1 R2>R3 yz=C2(1)+C2(2)*z+C2(3)*z.^2+C2(4)*z.^3+C2(5)*z.^5 "R2="; R2 endif if R3>R1 R3>R2 yz=C(1)+C(2)*z+C(3)*z.^3+C(4)*z.^5 "R3="; R3 endif % Построение графика. plot(x, y, "*r;experiment;",X1, Y1, ’-b;f(x);’,X1, Y2, ’dr;g(x);’,X1 , Y3, ’ok;fi(x);’, z, yz, ’sb;f(z);’); grid(); % Результаты работы программы. z = 1.0000 2.5000 4.8000 B = 0.083039 -0.567892 0.906779 1.609432 -1.115925 -1.355075 GGL = 1.3000e+01 2.8600e+01 1.5210e+02 7.2701e+02 1.2793e+05 2.8600e+01 1.5210e+02 7.2701e+02 3.9868e+03 7.5030e+05 1.5210e+02 7.2701e+02 3.9868e+03 2.2183e+04 4.4706e+06 7.2701e+02 3.9868e+03 2.2183e+04 1.2793e+05 2.6938e+07 2.2183e+04 1.2793e+05 7.5030e+05 4.4706e+06 1.6383e+08 A = 9.4262e+00 -3.6516e+00 -5.7767e+00 1.7888e+00 -5.8179e-05 C = -1.030345 5.080391 -0.609721 0.033534 R1 = 0.99690 R2 = 0.98136 R3 = 0.99573 yz = -0.43964 6.008544 0.77972
На рисунке 11.4 представлено графическое решение задачи.
(рис 11.4) Графическое решение к примеру 11.3
Рассмотренная задача демонстрирует основные приёмы подбора зависимости методом наименьших квадратов. Авторы рекомендует внимательно рассмотреть её для понимания методов решения подобных задач в Octave.
В заключении авторы позволят несколько советов по решению задачи аппроксимации.
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.