Функции Curve Fitting Toolbox

Простые примеры

Приближение без создания объекта с параметрической моделью

Предположим, что при помощи функций Curve Fitting Toolbox требуется приблизить некоторые данные полиномом третьей степени по методу наименьших квадратов. Сгенерируем эти данные и запишем их в вектор-столбцы (при созднии x используется ' для транспонирования и получения вектор-столбца из вектор-строки ). Массив данных содержит значения полинома третьей степени с ошибками, распределенными по нормальному закону со средним 0 и среднеквадратическим отклонением 1:

>> x=(0:0.05:3)';
>> y=5*x.^3-10*x.^2-3*x+randn(size(x));

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

• вектор-столбцы с данными;
• имя одной из моделей, в нашем случае для приближения полиномом третьей степени имя будет 'poly3' (стандартные параметрические и непараметрические модели описаны в разделах Параметрические модели библиотеки Curve Fitting Toolbox и Непараметрические модели библиотеки Curve Fitting Toolbox).

>> fresult = fit(x,y,'poly3')

Функция fit возвращает полученное приближение (в Curve Fitting Toolbox приближения и модели являются объектами).

fresult =
     Linear model Poly3:
       fresult(x) = p1*x^3 + p2*x^2 + p3*x + p4
     Coefficients (with 95% confidence bounds):
       p1 =       5.322  (4.916, 5.729)
       p2 =      -11.52  (-13.38, -9.667)
       p3 =      -1.064  (-3.447, 1.32)
       p4 =     -0.6663  (-1.485, 0.1524) 

В данном случае в командное окно MATLAB вывелся вид полиномиальной модели и полученные значения коэффициентов вместе с границами доверительных интервалов, соответствующих уровню вероятности 95%.

Для визуализации данных и построения графика полученного полинома третьей степени служит функция plot, переопределенная в Curve Fitting Toolbox. Если в качестве ее входного аргумента указать полученное приближение fresult, то построится его график. А если вызвать plot с тремя входными аргументами, указав еще и вектора с данными, то на графике отобразятся исходные данные:

>> plot(fresult, x, y)

График исходных данных и приближающего их полинома третьей степени

Переопределенная в Curve Fitting Toolbox функция plot удобна еще и тем, что имеет ряд дополнительных входных аргументов, при помощи которых можно построить, в частности: графики ошибок и доверительные полосы. Например следующее обращение к функции plot приводит к появлению в графическом окне еще и графика ошибок для построенного приближения:

>> plot(fresult,x,y,'residuals')

Графики данных, приближения и его ошибки

Если полученные значения коэффициентов параметрической модели требуется использовать в дальнейших расчетах, то их следует записать в некоторый массив, например coeff. Для этого нужно вызвать функцию coeffvalues , указав в качестве ее входного аргумента объект fresult:

>> coef=coeffvalues(fresult)
coef =
    4.9333   -9.3702   -4.2915    0.5525

Для вычисления значений параметрической модели в заданных точках, в принципе, можно пользоваться этими коэффициентами и функцией polyval, входящей в MATLAB. Однако, в составе Curve Fitting Toolbox имеется функция feval (это переопределенная функция feval, входящая в базовый набор функций MATLAB). Ей надо указать построенную параметрическую модель fresult и вектор с координатами точек, в которых требуется вычислить ее значения. Найдем, например, значения модели в точках от первого до последнего элемента вектора x с шагом 0.001:

>> xx=x(1):0.001:x(end);
>> ff = feval(fresult,xx);

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

>> plot(xx,ff) 

При таком варианте обращения к plot вызывается обычная функция plot, входящая в базовый набор функций MATLAB, а не переопределенная в Curve Fitting Toolbox.

Приближение с созданием объекта для параметрической модели

В приведенном выше примере мы сразу же вызвали функцию fit, входящую в Curve Fitting Toolbox, указав в ее входных аргументах исходные данные и имя параметрической модели, и получили приближение. Другой способ состоит в предварительном создании параметрической модели при помощи функции fittype и последующем нахождении ее параметров с использованием функции fit.

Приведем пример приближения данных:

>> x=(0:0.05:3)';
>> y=3.45*exp(-2.77*x)+0.05*randn(size(x));

при помощи экспоненциальной модели:

ae^bx

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

>> ftype = fittype('exp1')

В командное окно выводится вид этой модели:

ftype =
     General model Exp1:
       ftype(a,b,x) = a*exp(b*x)

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

>> fresult = fit(x,y,ftype)

В командное окно выводится информация об объекте fresult с полученным приближением: вид модели и найденные значения коэффициентов вместе с границами доверительных интервалов:

fresult =

     General model Exp1:
       fresult(x) = a*exp(b*x)
     Coefficients (with 95% confidence bounds):
       a =       3.471  (3.401, 3.542)
       b =      -2.754  (-2.838, -2.669)

Для визуализации данных и построения графика полученного приближения снова применим переопределенную в Curve Fitting Toolbox функцию plot:

plot(fresult,x,y)

получаем

Приближение данных моделью ae^bx

Приближение с созданием объекта для пользовательской линейной параметрической модели.

Функцию fittype можно применять и для создания собственной параметрической модели, в которую параметры входят как линейно, так и нелинейно. Рассмотрим сначала приближение данных при помощи линейной параметрической модели.

Сгенерируем исходные данные

>> x=(0:0.05:3)';
>> y=1.55*x+5*exp(-x)-2+0.05*randn(size(x));

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

p₁x + p₂e^-x + p₃

При создании объекта с линейной параметрической моделью при помощи функции fittype в качестве ее входного аргумента указывается массив ячеек (для создания массива ячеек используются фигурные скобки), каждая ячейка которого содержит строку со слагаемым, входящим в модель (без искомого коэффициента), т.е. {'x','exp(-x)','1'}. Если обратиться к функции fittype таким образом:

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

то получим линейную параметрическую модель, в которой коэффициентам автоматически даны имена a, b, c:

ftype =
     Linear model:
       ftype(a,b,c,x) = a*x + b*(exp(-x)) + c

Если же мы хотим, чтобы коэффициенты назывались именно p1, p2 и p3, то следует указать их имена при создании параметрической модели при помощи функции fittype. Для этого функция fittype вызывается с дополнительной парой аргументов: 'coefficients' и массивом ячеек, каждая ячейка которого содержит строку с именем соответствующего параметра модели {'p1','p2','p3'}:

>>  ftype = fittype({'x','exp(-x)','1'},'coefficients',{'p1','p2','p3'})

ftype =
     Linear model:
       ftype(p1,p2,p3,x) = p1*x + p2*(exp(-x)) + p3

Далее применяем функцию fit для подбора параметров

>> fresult = fit(x,y,ftype)

В командное окно выводится информация об объекте fresult с полученным приближением: вид модели и найденные значения параметров p₁, p₂ и p₃ вместе с границами доверительных интервалов, соответствующих уровню вероятности в 95%:

fresult =
     Linear model:
       fresult(x) = p1*x + p2*(exp(-x)) + p3
     Coefficients (with 95% confidence bounds):
       p1 =       1.543  (1.508, 1.578)
       p2 =       4.939  (4.821, 5.057)
       p3 =      -1.963  (-2.053, -1.873)

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

>> plot(fresult,x,y)

Приближение данных линейной моделью p₁x + p₂e^-x + p₃

Приближение с созданием объекта для пользовательской нелинейной параметрической модели.

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

>> x=(0:0.05:3)';
>> y=5.23*log(x+0.15)+0.5*randn(size(x));

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

aln(x + b)

Для этого создадим нелинейную параметрическую модель при помощи функции fittype, указав в ее входном аргументе строку 'a*log(x+b)' с выражением, задающим параметрическую модель:

>> ftype = fittype('a*log(x+b)')

Сама модель и сообщение о том, что она нелинейна, выводятся в командное окно:

ftype =
     General model:
       ftype(a,b,x) = a*log(x+b)

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

>> fresult = fit(x,y,ftype)

При этом выводится сообщение об ошибке, т.к. в процессе поиска параметров под знаком логарифма оказалось неположительное число.

??? Complex value computed by model function, fitting cannot continue.
Try using or tightening upper and lower bounds on coefficients. 

Для того, чтобы избежать таких ошибок, следует задавать границы допустимых интервалов изменения параметров. В нашем случае достаточно задать нижнюю границу только для второго параметра , скажем 1e-5. Для этого следует при помощи функции fitoptions сформировать управляющую структуру, которую затем необходимо указать во входных аргументах функции fit, выполняющей подбор параметров модели (более подробно об управлении вычислительными алгоритмами при подборе параметров написано в разд. Управление вычислительным алгоритмом). В функции fitoptions входные аргументы задаются парами:

НазваниеОпции, Значение, НазваниеОпции, Значение, …

Функция fitoptions возвращает управляющую структуру.

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

• Method - метод использующийся при нахождении приближения, в нашем случае это нелинейные наименьшие квадраты, т.е. значением опции Method должна быть строка 'NonlinearLeastSquares';
• Lower - вектор с нижними границами для значений параметров [-Inf 1e-5], где -Inf значит, что первый параметр может быть сколь угодно малым.

>> opts = fitoptions('Method','NonlinearLeastSquares','Lower',[-Inf 1e-5])

В командное окно выводится структура options

opts =
         Normalize: 'off'
          Exclude: []
          Weights: []
           Method: 'NonlinearLeastSquares'
           Robust: 'Off'
       StartPoint: [1x0 double]
            Lower: [-Inf 1.0000e-005]
            Upper: [1x0 double]
        Algorithm: 'Trust-Region'
    DiffMinChange: 1.0000e-008
    DiffMaxChange: 0.1000
          Display: 'Notify'
      MaxFunEvals: 600
          MaxIter: 400
           TolFun: 1.0000e-006
             TolX: 1.0000e-006

в которой поля Method и Lower имеют заданные значения. Смысл остальных полей объяснен в разделе Управление вычислительным алгоритмом.

Теперь укажем созданную управляющую структуру opts в четвертом входном аргументе функции fit

>> fresult = fit(x,y,ftype, opts)

В командное окно выводится информация об объекте fresult с полученным приближением: вид нелинейной параметрической модели и найденные значения параметров a и b вместе с границами доверительных интервалов, соответствующих уровню вероятности в 95%:

fresult =
     General model:
       fresult(x) = a*log(x+b)
     Coefficients (with 95% confidence bounds):
       a =       5.309  (5.141, 5.478)
       b =      0.1453  (0.1295, 0.1611)

Для визуализации результата осталось применить функцию plot:

>> plot(fresult,x,y)

Приближение данных пользовательской параметрической моделью aln(x + b)