Приближение данных, регрессия в ToolBox Optimization

Масштабирование данных (x_k,y_k)_k=1,2,...,n заключается замене x_k на , вычисляемые для k=1,2,...,n по формуле

где μ - среднее значение

а σ соответственно - среднеквадратичное отклонение

Для автоматического  масштабирования данных требуется  при формировании функцией fitoptions объекта opts с опциями вычислительного алгоритма указать для опции Normalize значение 'on' (по умолчанию опция Normalize принимает значение 'off', т.е. масштабирование и центрирование не применяется).

Приведем пример. Предположим, что требуется найти приближение следующего набора данных:

>> xdata = [1951;  1952;  1953;  1954;  1955;  1956;  1957];

>> ydata = [1.2; 3.4; 2.9; 4.4; 4.5; 5.1; 4.2];

полиномом четвертой  степени. Имя такой полиномиальной модели poly4 (стандартные параметрические  и непараметрические модели, входящие в библиотеку моделей Curve Fitting Toolbox описаны в разд. Параметрические и непараметрические модели). Создадим соответствующий ей объект ftype при помощи функции fittype:

>> ftype=fittype('poly4');

Приблизим сначала  данные при помощи функции fit без  предварительного центрирования и  масштабирования (как делается по умолчанию), найденное приближение запишем  в объект fresult:

>> fresult=fit(xdata,ydata,ftype);

При этом в командное  окно выводится сообщение о том, что соответствующая система  линейных уравнений плохообусловлена:

Warning: Equation is badly conditioned. Remove repeated data points

         or try centering and scaling.

Выведем найденные  значения параметров модели (т.е. коэффициенты полинома четвертой степени) при помощи функции coeffvalues, входящей в Curve Fitting Toolbox. Ее входным аргументом является объект fresult с найденным приближением:

>> format short e

>> c=coeffvalues(fresult)

 

c =

-3.5997e-002  2.8134e+002 -8.2459e+005  1.0741e+009 -5.2470e+011

Посмотрим теперь, как изменятся параметры модели, если в исходные данные внести небольшую  ошибку. Добавим к y(1) число 0.01 и измененный вектор значений запишем в переменную ydataerr

>> ydataerr = [1.21; 3.4; 2.9; 4.4; 4.5; 5.1; 4.2];

Снова выполним приближение данных xdata, ydataerr при  помощи полинома четвертой степени  без предварительного центрирования  и масштабирования и выведем  получающиеся значения параметров:

>> fresult=fit(xdata,ydataerr,ftype)

>> cerr=coeffvalues(fresult)

 

cerr =

-3.5883e-002  2.8045e+002 -8.2198e+005  1.0707e+009 -5.2304e+011

Видим, что ошибка в исходных данных в 0.01 привела к  ошибке порядка 1e9 в коэффициентах  полинома.

Теперь проделаем  то же самое с предварительным  центрированием и масштабированием данных, сформировав при помощи функции fittype объект opts со включенной опцией центрирования и масштабирования:

>> opts=fitoptions('Normalize','on');

>> nmfresult = fit(xdata,ydata,ftype,opts);

>> c=coeffvalues(nmfresult)

>> nmfresult = fit(xdata,ydataerr,ftype,opts);

>> cerr=coeffvalues(nmfresult)

 

c =

-7.8367e-001 -8.4010e-002  8.9621e-001  1.2061e+000  3.9108e+000

 

cerr =

-7.8120e-001 -8.6810e-002  8.9391e-001  1.2080e+000  3.9110e+000

Итак, предварительное  центрирование и масштабирование данных позволило получить ошибку в параметрах того же порядка, что и ошибка в исходных данных.

Опции для интерполяционных сплайнов

Если в Curve Fitting Toolbox производится интерполяция:

• кусочно-постоянными функциями по ближайшему соседу (опции method установлено значение NearestInterpolant);
• кусочно-линейными функциями (опции method установлено значение LinearInterpolant);
• кубическими сплайнами (опции method установлено значение CubicSplineInterpolant);
• сплайнами, сохраняющими форму данных, т.е. кусочными полиномами Эрмита (опции method установлено значение PchipInterpolant);

то никаких  дополнительных опций задавать нельзя и к функции fitoptions обращаются, соответственно, одним из четырех следующих способов:

opts=fitoptions('method', 'NearestInterpolant')

opts=fitoptions('method', 'LinearInterpolant')

opts=fitoptions('method', 'CubicSplineInterpolant')

opts=fitoptions('method', 'PchipInterpolant')

Такой вариант  задания опций функцией fitoptions практически  не используется, поскольку перечисленные способы интерполяции могут быть указаны сразу при создании объекта с моделью, например, для интерполяции кубическими полиномами Эрмита табличной функции:

>> xdata=(0:0.3:3)';

>> ydata=sin(xdata).^2.*exp(-xdata);

достаточно  выполнить следующие команды:

 

>> ftype=fittype('pchipinterp');

>> fresult=fit(xdata,ydata,ftype)

>> plot(fresult,xdata,ydata)

 
Интерполяция кубическими полиномами Эрмита

Опции для сглаживающего сплайна

Если требуется  приблизить данные сглаживающим сплайном, то опции method должно быть установлено  значение SmoothingSpline. Тогда появляется возможность задать значение параметра, входящего в сглаживающий сплайн. Для этого надо установить опции SmoothingParam необходимое значение от 0 до 1.

Например, приблизим  данные

>> xdata=(0:0.05:3)';

>> ydata=xdata.*exp(-xdata)+0.01*randn(size(xdata));

сглаживающим  сплайном со значением его параметра 0.9. Для этого сформируем объект opts, задав соответствующие опции во входных аргументах функции fitoptions:

>> opts = fitoptions('Method','SmoothingSpline','SmoothingParam',0.9)

Затем сформируем объект ftype с моделью при помощи функции fittype, указав название модели в  ее входном аргументе (стандартные параметрические и непараметрические модели Curve Fitting Toolbox описаны в разделе Параметрические и непараметрические модели)

>> ftype = fittype('smoothingspline')

Вызовем функцию fit для получения приближения fresult и применим переопределенную в Curve Fitting Toolbox функцию plot для визуализации исходных данных и вывода графика полученного сглаживающего сплайна

>> fresult = fit(xdata,ydata,ftype,opts)

>> plot(fresult,xdata,ydata)

 
Приближение данных сглаживающим сплайном

Опции для линейных параметрических моделей

Примечание.

Описанные в  этом разделе опции задаются также  и для нелинейных параметрических моделей.

В случае, когда  параметрическая модель содержит только линейно входящие параметры, набор  возможных опций, управляющих вычислительным алгоритмом, состоит из выбора устойчивого  по отношению к выбросам в данных алгоритма нахождения параметров модели и указания верхних и нижних границ для искомых параметров линейной параметрической модели:

• Robust - использование устойчивого процесса поиска параметров модели по отношению к выбросам в данных, возможные значения этой опции: 'on' (использовать устойчивый метод), 'off' (по умолчанию устойчивый метод не используется).
• Lower - Вектор из значений нижних границ для искомых параметров линейной модели. Число элементов в нем должно совпадать с числом параметров модели, если некоторые параметры могут быть сколь угодно малы, то на соответствующей позиции в векторе ставится -Inf (в MATLAB так обозначается минус бесконечность). Если ограничения снизу на параметры модели отсутствуют, то значением является пустой массив, т.е [] (по умолчанию).
• Upper - Вектор из значений верхних границ для искомых параметров линейной модели. Число элементов в нем должно совпадать с числом параметров модели, если некоторые параметры могут быть сколь угодно велики, то на соответствующей позиции в векторе ставится Inf (в MATLAB так обозначается плюс бесконечность). Если ограничения сверху на параметры модели отсутствуют, то значением является пустой массив, т.е [] (по умолчанию).

В векторах, являющихся значениями опций Lower и Upper, ограничения  должны задаваться так, чтобы они соответствовали упорядоченным по алфавиту коэффициентам параметрической модели, т.е., если линейная пользовательская параметрическая модель имеет вид:

то, например вектор [0.1 2 -3 4.8], являющийся значением свойства Lower, соответствует следующим ограничениям снизу на параметры модели:

Наличие выбросов в данных может сильно ухудшить качество приближения параметрической моделью, поскольку при поиске параметров исходя из условия минимизации суммы квадратов ошибок

где (x_k,y_k)_k=1,2,...,n - исходные данные, расстояние от искомой кривой до выброса окажет сильное влияние на значения параметров модели.

Одной из возможностей автоматического снижения влияния  выбросов является назначение маленького веса соответствующим данным, или полное исключение выбросов из процесса подбора параметров (см. разделы Опция для задания весов и Опция для исключения части данных). Однако, такая ручная предварительная обработка не всегда может привести к желаемому результату, поскольку не всегда можно точно идентифицировать выброс, или подобрать нужный вес для снижения его влияния на приближение. Поэтому в Curve Fitting Toolbox предусмотрен адаптивный алгоритм построения приближения, в ходе которого данным автоматически присваиваются веса в зависимости от их удаления от полученной на данном шаге линии графика параметрической модели. Для задействования адаптивного алгоритма с целью снижения влияния выбросов в данных на приближение следует установить опции Robust значение 'on'.

В адаптивном алгоритме  для снижения влияния выбросов на каждом шаге строится приближение по методу наименьших квадратов и вычисляются приведенные невязки u_i (см. более подробно в справочной системе по Curve Fitting Toolbox раздел Fitting Data: Parametric Fitting: The Least Squares Fitting Method) и в зависимости от удаленности точек от приближающей их кривой данным назначаются веса по следующему правилу (чем дальше точка, тем меньше ее вес):

Далее процесс  приближения повторяется, пока не выполнится некоторое условие сходимости.

Поскольку в  описанном выше алгоритме при  вычислении весов используется четвертая  степень, то он называется bisquare weights (биквадратные веса).

Приведем пример, в котором задание опции Robust значения 'on' улучшает качество приближения. Рассмотрим следующие исходные данные, полученные из закона

y = 3.8x + 2.1

внесением распределенной по нормальному закону ошибки

>> xdata=(0:0.5:11)';

>> ydata=3.8*xdata+2.1+2*randn(size(xdata));

в которых есть два выброса (см. график ниже)

>> ydata(2)=35;

>> ydata(20)=5;

 
Данные с выбросами

При помощи функции fittype построим линейную параметрическую  модель - полином первой степени - для проведения линейной регрессии (создание пользовательских линейных моделей в Curve Fitting Toolbox описано в разделе Приближение с созданием объекта для пользовательской линейной параметрической модели).

>> ftype=fittype({'x','1'})

 

ftype =

     Linear model:

       ftype(a,b,x) = a*x + b

Приблизим сначала  данные при помощи обычного метода наименьших квадратов, не указывая никаких  опций (тогда обычный метод наименьших квадратов применяется по умолчанию):

>> fres=fit(xdata,ydata,ftype)

Получаем следующие  значения параметров модели

res =

     Linear model:

       fres(x) = a*x + b

     Coefficients (with 95% confidence bounds):

       a =       2.638  (1.421, 3.854)

       b =       8.297  (0.4845, 16.11)

Теперь создадим объект с нужными настройками  вычислительного алгоритма:

• опция Method должна принимать значение 'LinearLeastSquares';
• опция Robust должна принимать значение 'on'.

>> opts=fitoptions('Method', 'LinearLeastSquares','Robust','on');

Приблизим данные с выбросами при помощи адаптивного  алгоритма для снижения влияния  выбросов на получающееся приближение:

>> frobres=fit(xdata,ydata,ftype,opts)

 

frobres =

     Linear model:

       frobres(x) = a*x + b

     Coefficients (with 95% confidence bounds):

       a =       3.739  (3.199, 4.28)

       b =       2.262  (-1.21, 5.734)

Полученные  при помощи адаптивного алгоритма  результаты соответствуют исходным данным, в отличие от тех, что были получены обычным методом наименьших квадратов, в котором по умолчанию всем данным приписаны равные единице веса. Визуальное сравнение двух приближений также позволяет в этом убедиться. Построим в одном графическом окне на разных осях (создаваемых при помощи функции subplot) приближения вместе с исходными данными

>> figure

>> subplot(1,2,1)

>> plot(fres,xdata,ydata)

>> subplot(1,2,2)

>> plot(frobres,xdata,ydata)

В результате получаем