В этом примере показаны подводные камни, которые могут возникнуть при подгонке нелинейной модели путем преобразования в линейность. Представьте, что мы собрали измерения по двум переменным, 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);

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