Стратегия подгонки сильно нелинейной функции

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

$y = ax + bx^{-1/2}$

Здесь, особенно значение представляет большой интерес. $b$

Сюжет для этой функции:

Функциональный сюжет

(Обратите внимание, что функция модели основана на подробном математическом описании системы и, кажется, работает очень хорошо - просто автоматические подборы сложны).

Конечно, функция модели проблематична: подходящие стратегии, которые я пробовал до сих пор, терпят неудачу из-за резкой асимптоты при , особенно с зашумленными данными. $x=0$

Мое понимание проблемы здесь заключается в том, что простая аппроксимация методом наименьших квадратов (я играл как с линейной, так и с нелинейной регрессией в MATLAB; в основном это Левенберг-Марквардт) очень чувствительна к вертикальной асимптоте, потому что небольшие ошибки в x сильно усиливаются ,

Может ли кто-нибудь указать мне подходящую стратегию, которая может обойти это?

У меня есть некоторые базовые знания в области статистики, но они все еще довольно ограничены. Я был бы готов учиться, если бы я только знал, где начать искать :)

Большое спасибо за ваш совет!

Править Прошу прощения за то, что забыл упомянуть об ошибках. Единственный значительный шум в , и он аддитивен. $x$

Изменить 2 Некоторые дополнительные сведения об истории вопроса. График выше моделирует поведение при растяжении полимера. Как отметил @whuber в комментариях, вам нужно чтобы получить график, как показано выше. $b \approx -200 a$

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

Редактировать 3 и 4 Фиксированный график.

curve-fitting nonlinear onnodb
источник

Есть ошибки в или или обоих? В какой форме вы ожидаете шум (мультипликативный, аддитивный и т. Д.)?

x

$x$

y

$y$

вероятностная

@onnodb: Меня беспокоит, разве это не принципиальный вопрос, насколько надежна ваша модель? Независимо от того, какую подходящую стратегию вы используете, не останется очень чувствительным? Можете ли вы когда-нибудь быть уверенным в такой оценке для ?

b

$b$

b

$b$

curious_cat

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

a

$a$

b

$b$

b

$b$

a

$a$

a

$a$

b x^{1 / 2}

$b x^{1/2}$

whuber

Спасибо, но это все еще не так. Расширяя касательную к этому графику в обратном направлении из любой точки где , вы ось y в точке . Поскольку нисходящий всплеск в показывает, что отрицателен, этот y-перехват также должен быть отрицательным. Но на вашей фигуре совершенно очевидно, что большинство таких перехватов являются положительными, вплоть до . Таким образом, математически невозможно, чтобы уравнение типа могло описать вашу кривую , даже приблизительно. Как минимум, вам нужно подобрать что-то вроде .

(x, a x + b x^{1 / 2})

$(x,ax+bx^{1/2})$

x > 0

$x\gt 0$

(0, 3 b / (2 x^{1 / 2}))

$(0,3b/(2x^{1/2}))$

0

$0$

b

$b$

15.5

$15.5$ $y=ax + bx^{1/2}$

y = a x + b x^{1 / 2} + c

$y=ax+bx^{1/2}+c$

whuber

Прежде чем я поработал над этим, я хотел убедиться в правильности постановки вопроса: вот почему важно правильно настроить функцию. Сейчас у меня нет времени, чтобы дать полный ответ, но я хотел бы отметить, что «другие люди» могут ошибаться, но, увы, это зависит от еще большего количества деталей. Если ваша ошибка действительно аддитивна, мне кажется, она все равно должна быть сильно гетероскедастичной, иначе ее дисперсия при малых значениях будет действительно крошечной. Что вы можете сказать нам - количественно - об этой ошибке?

x

$x$

x

$x$

whuber

Ответы:

Методы, которые мы использовали бы для подбора этого вручную (то есть исследовательского анализа данных), могут замечательно работать с такими данными.

Я хочу немного изменить параметры модели , чтобы сделать ее параметры положительными:

y = a x - b / \sqrt{x} .

$y = a x - b / \sqrt{x}.$

Для данного , давайте предположим, что существует единственное действительное удовлетворяющее этому уравнению; назовите это или, для краткости, когда поняты. $y$ $x$ $f(y; a,b)$ $f(y)$ $(a,b)$

Мы наблюдаем набор упорядоченных пар где отклоняются от независимыми случайными переменными с нулевым средним. В этом обсуждении я предполагаю, что все они имеют общую дисперсию, но расширение этих результатов (с использованием взвешенных наименьших квадратов) возможно, очевидно и легко реализуемо. Вот смоделированный пример такого набора из значений с , и общей дисперсией . $(x_i, y_i)$ $x_i$ $f(y_i; a,b)$ $100$ $a=0.0001$ $b=0.1$ $\sigma^2=4$

График данных

Это (намеренно) жесткий пример, который можно оценить по нефизическим (отрицательным) значениям и их необычайному разбросу (который обычно равен горизонтальным единицам, но может варьироваться до или по оси ). Если мы сможем получить разумное соответствие этим данным, которое приблизится к оценке используемых , и , мы действительно преуспеем. $x$ $\pm 2$ $5$ $6$ $x$ $a$ $b$ $\sigma^2$

Поисковая примерка носит итеративный характер. Каждая стадия состоит из двух этапов: оценки (на основе данных , и предыдущих оценки и из и , из которого предыдущих предсказанных значений может быть получены для ), а затем оценить . Поскольку ошибки в x , подгонки оценивают по , а не наоборот. К первому порядку ошибок по , когда достаточно большой, $a$ $\hat{a}$ $\hat{b}$ $a$ $b$ $\hat{x}_i$ $x_i$ $b$ $x_i$ $(y_i)$ $x$ $x$

x_{i} \approx \frac{1}{a} (y_{i} + \frac{\hat{b}}{\sqrt{{\hat{x}}_{i}}}) .

$x_i \approx \frac{1}{a}\left(y_i + \frac{\hat{b}}{\sqrt{\hat{x}_i}}\right).$

Поэтому мы можем обновить , подгоняя эту модель наименьшими квадратами (обратите внимание, что она имеет только один параметр - наклон, - и не перехватывая) и принимая обратную величину коэффициента в качестве обновленной оценки . $\hat{a}$ $a$ $a$

Затем, когда достаточно мало, обратный квадратичный член доминирует, и мы находим (снова в первом порядке по ошибкам), что $x$

x_{i} \approx b^{2} \frac{1 - 2 \hat{a} \hat{b} {\hat{x}}^{3 / 2}}{y_{i}^{2}} .

$x_i \approx b^2\frac{1 - 2 \hat{a} \hat{b} \hat{x}^{3/2}}{y_i^2}.$

Еще раз, используя метод наименьших квадратов (только с углом наклона ), мы получаем обновленную оценку через корень квадратный из подогнанного наклона. $b$ $\hat{b}$

Чтобы понять, почему это работает, грубая исследовательская аппроксимация этой подгонки может быть получена путем построения графика против для меньшего . Еще лучше, потому что измеряются с ошибкой, а монотонно изменяется с , мы должны сосредоточиться на данных с большими значениями . Вот пример из нашего смоделированного набора данных, показывающий наибольшую половину красным, наименьшую половину синим и линию через начало координат, подходящую к красным точкам. $x_i$ $1/y_i^2$ $x_i$ $x_i$ $y_i$ $x_i$ $1/y_i^2$ $y_i$

фигура

Точки приблизительно совпадают, хотя при малых значениях и наблюдается небольшая кривизна . (Обратите внимание на выбор осей: поскольку является измерением, условно нанести его на вертикальную ось.) Сосредоточив посадку на красных точках, где кривизна должна быть минимальной, мы должны получить разумную оценку . Значение указанное в заголовке, является квадратным корнем наклона этой линии: это только на % меньше истинного значения! $x$ $y$ $x$ $b$ $0.096$ $4$

На данный момент прогнозируемые значения могут быть обновлены с помощью

{\hat{x}}_{i} = f (y_{i}; \hat{a}, \hat{b}) .

$\hat{x}_i = f(y_i; \hat{a}, \hat{b}).$

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

Оказывается, что трудно оценить, если у нас нет хорошего набора очень больших значений , но который определяет вертикальную асимптоту в исходном графике (в вопросе) и является предметом вопроса - может быть точно определено, если в вертикальной асимптоте есть некоторые данные. В нашем рабочем примере итерации сходятся к (что почти вдвое больше правильного значения ) и (что близко к правильному значению ). Этот график еще раз показывает данные, на которые накладываются (а) истинные $a$ $x$ $b$ $\hat{a} = 0.000196$ $0.0001$ $\hat{b} = 0.1073$ $0.1$ кривая серым цветом (пунктирная) и (b) расчетная кривая красным цветом (сплошная):

Приступы

Эта подгонка настолько хороша, что трудно отличить истинную кривую от подгоночной кривой: они перекрываются почти везде. Кстати, оценочная дисперсия ошибки очень близка к истинному значению . $3.73$ $4$

Есть некоторые проблемы с этим подходом:

Оценки предвзяты. Смещение становится очевидным, когда набор данных небольшой, и относительно немного значений находятся близко к оси x. Подгонка систематически немного низкая.
Процедура оценки требует, чтобы метод отличал "большие" от "маленьких" значений . Я мог бы предложить исследовательские способы определения оптимальных определений, но на практике вы можете оставить их как «настраивающие» константы и изменить их, чтобы проверить чувствительность результатов. Я установил их произвольно, разделив данные на три равные группы в соответствии со значением и используя две внешние группы. $y_i$ $y_i$
Процедура не будет работать для всех возможных комбинаций и или всех возможных диапазонов данных. Однако он должен работать хорошо, когда в наборе данных представлено достаточно кривой, чтобы отразить обе асимптоты: вертикальную на одном конце и наклонную на другом конце. $a$ $b$

Код

Следующее написано в Mathematica .

estimate[{a_, b_, xHat_}, {x_, y_}] := 
  Module[{n = Length[x], k0, k1, yLarge, xLarge, xHatLarge, ySmall, 
    xSmall, xHatSmall, a1, b1, xHat1, u, fr},
   fr[y_, {a_, b_}] := Root[-b^2 + y^2 #1 - 2 a y #1^2 + a^2 #1^3 &, 1];
   k0 = Floor[1 n/3]; k1 = Ceiling[2 n/3];(* The tuning constants *)
   yLarge = y[[k1 + 1 ;;]]; xLarge = x[[k1 + 1 ;;]]; xHatLarge = xHat[[k1 + 1 ;;]];
   ySmall = y[[;; k0]]; xSmall = x[[;; k0]]; xHatSmall = xHat[[;; k0]];
   a1 = 1/
     Last[LinearModelFit[{yLarge + b/Sqrt[xHatLarge], 
          xLarge}\[Transpose], u, u]["BestFitParameters"]];
   b1 = Sqrt[
     Last[LinearModelFit[{(1 - 2 a1 b  xHatSmall^(3/2)) / ySmall^2, 
          xSmall}\[Transpose], u, u]["BestFitParameters"]]];
   xHat1 = fr[#, {a1, b1}] & /@ y;
   {a1, b1, xHat1}
   ];

Примените это к данным (заданным параллельными векторами xи yсформированными в матрицу из двух столбцов data = {x,y}) до сходимости, начиная с оценок : $a=b=0$

{a, b, xHat} = NestWhile[estimate[##, data] &, {0, 0, data[[1]]}, 
                Norm[Most[#1] - Most[#2]] >= 0.001 &,  2, 100]

Whuber
источник

Это удивительный ответ; Я очень благодарен! Я играл с этим, и результаты выглядят очень многообещающе. Хотя мне понадобится немного больше времени, чтобы полностью понять причину :) Кроме того: могу ли я связаться с вами через ваш веб-сайт для еще одного (частного) вопроса относительно подтверждений?

onnodb

Смотрите важные вопросы @probabilityislogic опубликовал

Если у вас есть только ошибки в y, и они аддитивны, и у вас постоянная дисперсия (т.е. ваши предположения соответствуют тому, что звучит так, как вы это делали), то если вы позволите , вы, возможно, попытаетесь взвешенная линейная аппроксимация на , где веса будут пропорциональны ... (и да, это может просто сдвинуть проблему, так что может все еще быть проблематичным - но вы должны, по крайней мере, легче упорядочить с этим преобразованием проблемы). $y^* = y\sqrt{x}$ $y^*$ $x^* = x^{3/2}$ $1/x$

Обратите внимание, что с этой манипуляцией ваш становится перехватом нового уравнения $b$

Если ваши отклонения уже не постоянны, или ваши ошибки не являются аддитивными, или у вас есть ошибки в , это изменит вещи. $x$

Изменить, чтобы рассмотреть дополнительную информацию:

Мы добрались до модели вида: $y^* = b + a x^*$

Теперь у нас есть ошибки в х и аддитивные. Мы до сих пор не знаем, является ли дисперсия постоянной в этом масштабе.

Перепишите как $x^* = y^*/a - b/a = m y^* + c$

Пусть , где этот термин ошибки может быть гетероскедастичным (если исходный имеет постоянный разброс, он будет гетероскедастичным, но известной формы) $x_o^* = x^* + \eta$ $x$

(где в означает «наблюдаемый») $o$ $x_o^*$

Тогда где выглядит хорошо, но теперь имеет коррелированные ошибки в переменных и ; так что это линейная модель ошибок в переменных, с гетероскедастичностью и известной формой зависимости в ошибках. $x^*_o = c + m y^* + \epsilon$ $\epsilon = -\zeta$ $x$ $y$

Я не уверен, что это улучшает ситуацию! Я считаю, что есть методы для такого рода вещей, но это не совсем моя сфера.

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

Возможно, вы даже застряли с попыткой использовать методы, устойчивые к ошибкам в этой линейной форме.

Теперь огромный вопрос: если ошибки в х , как, черт возьми, вы подходите нелинейной модели? Вы просто слепо минимизировали сумму квадратов ошибок по ? Это может быть вашей проблемой. $y$

Я полагаю, что можно попытаться переписать исходную вещь как модель с ошибками в и попытаться оптимизировать подгонку, но я не уверен, что вижу, как правильно это настроить. $x$

Glen_b - Восстановить Монику
источник

Благодарность! Это интересное преобразование, я об этом даже не думал - хотя ошибки в , я все равно поиграю с этим!

x

$x$

onnodb

« Даже если ошибки в х » - это очень важно. Вы можете проверить обратную регрессию.

Glen_b

... или вы могли бы соответствовать модели :-).

x = \frac{1}{3} (\frac{2 y}{a} + \frac{2^{1 / 3} y^{2}}{{(27 a^{4} b^{2} - 2 a^{3} y^{3} + 3 \sqrt{3} \sqrt{27 a^{8} b^{4} - 4 a^{7} b^{2} y^{3}})}^{1 / 3}} + \frac{{(27 a^{4} b^{2} - 2 a^{3} y^{3} + 3 \sqrt{3} \sqrt{27 a^{8} b^{4} - 4 a^{7} b^{2} y^{3}})}^{1 / 3}}{2^{1 / 3} a^{2}})

$x = \frac{1}{3} \left(\frac{2 y}{a}+\frac{2^{1/3} y^2}{\left(27 a^4 b^2-2 a^3 y^3+3 \sqrt{3} \sqrt{27 a^8 b^4-4 a^7 b^2 y^3}\right)^{1/3}}+\frac{\left(27 a^4 b^2-2 a^3 y^3+3 \sqrt{3} \sqrt{27 a^8 b^4-4 a^7 b^2 y^3}\right)^{1/3}}{2^{1/3} a^2}\right)$

whuber

@ Whuber Хмм. Решаю за куб, умный. Если мы напишем оригинал в терминах где - , это оставит нас с (снова с ), что хотя бы условно можно сделать с помощью нелинейных наименьших квадратов. Похоже, что он позаботился о распространении ошибок должным образом. Это может на самом деле сработать, если OP будет использовать линейную форму, с которой я играл (используя некоторую робастную оценку ошибок-в-IV-и-гетеро), чтобы получить хорошие начальные значения для параметров, а затем попытаться использовать эта нелинейная форма LS, чтобы полировать это.

x_{o}

$x_o$

x_{o}

$x_o$

x + ζ

$x + \zeta$

x = (t h a t m o n s t e r) + ϵ

$x = (\rm{that\,\, monster}) + \epsilon\,$

ϵ = - ζ

$\epsilon = -\zeta$

Glen_b

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

x (y)

$x(y)$

y

$y$

b

$b$

whuber

После еще нескольких недель экспериментов, в данном конкретном случае лучше всего работает другая техника: подгонка Total Least Squares . Это вариант обычной (нелинейной) подгонки наименьших квадратов, но вместо измерения ошибок подгонки только по одной из осей (что вызывает проблемы в сильно нелинейных случаях, таких как эта), он учитывает обе оси.

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

onnodb
источник

Спасибо, что поделились этим. Я согласен, что это может привести к хорошим результатам в вашем случае, но у меня есть две проблемы. Первое, что вы упомянули: как точно применить метод наименьших квадратов / регрессию ошибок в переменных / ортогональную регрессию / регрессию Деминга к нелинейным подгонкам? Во-вторых, этот подход не подходит для ваших данных, в которых измеряется практически без ошибок. В этом случае вы не должны допускать наличие остатков в переменной и это должно привести к ненадежным, смещенным результатам.

y

$y$

y

$y$

whuber

@whuber Спасибо за выражение ваших проблем! Прямо сейчас я все еще работаю над запуском симуляций, чтобы исследовать надежность подгонки TLS для этой проблемы. Однако до сих пор я видел, что рассмотрение обеих переменных в TLS очень помогает в преодолении высокой нелинейности модели. Подходы смоделированных данных надежны и сходятся очень хорошо. Однако необходимо проделать еще больше работы, и мне непременно придется довести ваш метод до этого, как только у нас появится больше фактических данных, - и подробно рассмотреть ваши проблемы.

onnodb

Хорошо - не забывайте, у меня есть сопоставимые проблемы относительно метода, который я предложил!

whuber