Документация

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