Этот пример показывает подводные камни, которые могут возникнуть при подборе нелинейной модели путем преобразования в линейность. Представьте, что мы собрали измерения по двум переменным, x и y, и мы хотим смоделировать y как функцию x. Предположим, что x измеряется точно, в то время как измерения y затронуты аддитивными, симметричными, нулевыми ошибками.

x = [5.72 4.22 5.72 3.59 5.04 2.66 5.02 3.11 0.13 2.26 ...
     5.39 2.57 1.20 1.82 3.23 5.46 3.15 1.84 0.21 4.29 ...
     4.61 0.36 3.76 1.59 1.87 3.14 2.45 5.36 3.44 3.41]';
y = [2.66 2.91 0.94 4.28 1.76 4.08 1.11 4.33 8.94 5.25 ...
     0.02 3.88 6.43 4.08 4.90 1.33 3.63 5.49 7.23 0.88 ...
     3.08 8.12 1.22 4.24 6.21 5.48 4.89 2.30 4.13 2.17]';

Давайте также предположим, что теория говорит нам, что эти данные должны следовать модели экспоненциального распада, y = p1 * exp (p2 * x), где p1 положительно, а p2 отрицательно. Чтобы соответствовать этой модели, мы могли бы использовать нелинейные методы наименьших квадратов.

modelFun = @(p,x) p(1)*exp(p(2)*x);

Но нелинейная модель также может быть преобразована в линейную путем взятия логарифмической с обеих сторон, чтобы получить логарифмическую (y) = логарифмическую (p1) + p2 * x. Это заманчиво, потому что мы можем подгонять эту линейную модель обычными линейными методами наименьших квадратов. Коэффициенты, которые мы получили бы из линейного метода наименьших квадратов, будут журналом (p1) и p2.

paramEstsLin = [ones(size(x)), x] \ log(y);
paramEstsLin(1) = exp(paramEstsLin(1))

paramEstsLin =

   11.9312
   -0.4462

Как мы это сделали? Мы можем наложить подгонку на данные, чтобы выяснить.

xx = linspace(min(x), max(x));
yyLin = modelFun(paramEstsLin, xx);
plot(x,y,'o', xx,yyLin,'-');
xlabel('x'); ylabel('y');
legend({'Raw data','Linear fit on the log scale'},'location','NE');

Похоже, что-то пошло не так, потому что подгонка на самом деле не соответствует тренду, который мы видим в необработанных данных. Какую подгонку мы получили бы, если бы использовали nlinfit вместо этого делать нелинейные методы наименьших квадратов? Мы будем использовать предыдущую подгонку как грубую начальную точку, хотя это не очень подходит.

paramEsts = nlinfit(x, y, modelFun, paramEstsLin)

paramEsts =

    8.8145
   -0.2885

yy = modelFun(paramEsts,xx);
plot(x,y,'o', xx,yyLin,'-', xx,yy,'-');
xlabel('x'); ylabel('y');
legend({'Raw data','Linear fit on the log scale',  ...
	'Nonlinear fit on the original scale'},'location','NE');

Подгонка с использованием nlinfit более или менее проходит через центр рассеяния точек данных. Остаточный график показывает что-то примерно похожее на ровное рассеяние около нуля.

r = y-modelFun(paramEsts,x);
plot(x,r,'+', [min(x) max(x)],[0 0],'k:');
xlabel('x'); ylabel('residuals');

Так что же пошло не так с линейной подгонкой? Задача заключается в преобразовании журнала. Если мы постройте график данных и двух подгонка по шкале журнала, мы увидим, что есть экстремальные выбросы.

plot(x,log(y),'o', xx,log(yyLin),'-', xx,log(yy),'-');
xlabel('x'); ylabel('log(y)');
ylim([-5,3]);
legend({'Raw data', 'Linear fit on the log scale',  ...
	'Nonlinear fit on the original scale'},'location','SW');

Это наблюдение не является выбросами в исходных данных, так что же случилось, чтобы сделать его по журналу шкале? Преобразование журнала - это именно то, что нужно для выпрямления линии тренда. Но журнал является очень нелинейным преобразованием, и поэтому симметричные ошибки измерения на исходной шкале стали асимметричными на журнал шкале. Заметьте, что выбросы имели наименьшее значение y по исходной шкале - близкое к нулю. Преобразование журнала «растянуло», что наименьшее значение y больше, чем его соседи. Мы сделали линейную подгонку на журнал шкале, и поэтому это очень сильно затронуто этими выбросами.

Если бы измерение в этой одной точке было немного иным, то две подгонки могли бы быть намного более похожими. Для примера,

y(11) = 1;
paramEsts = nlinfit(x, y, modelFun, [10;-.3])

paramEsts =

    8.7618
   -0.2833

paramEstsLin = [ones(size(x)), x] \ log(y);
paramEstsLin(1) = exp(paramEstsLin(1))

paramEstsLin =

    9.6357
   -0.3394

yy = modelFun(paramEsts,xx);
yyLin = modelFun(paramEstsLin, xx);
plot(x,y,'o', xx,yyLin,'-', xx,yy,'-');
xlabel('x'); ylabel('y');
legend({'Raw data', 'Linear fit on the log scale',  ...
	'Nonlinear fit on the original scale'},'location','NE');

Тем не менее, этот два подгонка разные. Кто из них «прав»? Чтобы ответить на это, предположим, что вместо аддитивных ошибок измерения y были затронуты мультипликативными ошибками. Эти ошибки не будут симметричными, и наименьшие квадраты в исходной шкале не будут подходящими. С другой стороны, преобразование журнала сделает ошибки симметричными на шкале журналами, и линейном методе наименьших квадратов, подгонке на этой шкале, подходит.

Итак, какой метод «прав», зависит от того, какие предположения вы готовы сделать о своих данных. На практике, когда термин шума является маленьким относительно тренда, журнал преобразование является «локально линейным» в том смысле, что значения y вблизи того же значения x не будут растянуты слишком асимметрично. В этом случае эти два метода приводят по существу к одинаковой подгонке. Но когда термин шума не является маленьким, следует учитывать, какие предположения являются реалистичными, и выбрать соответствующий метод подгонки.