擬合高度非線性函數的策略

March 13, 2013

為了分析來自生物物理學實驗的數據，我目前正在嘗試使用高度非線性模型進行曲線擬合。模型函數看起來基本上像：

在這裡，尤其是價值非常感興趣。

此函數的繪圖：

（請注意，模型函數基於對系統的全面數學描述，並且似乎工作得很好——只是自動擬合很棘手）。

當然，模型函數是有問題的：到目前為止我嘗試過的擬合策略，由於尖銳的漸近線而失敗，尤其是有噪聲的數據。

我對這裡問題的理解是簡單的最小二乘擬合（我在 MATLAB 中使用過線性和非線性回歸；主要是 Levenberg-Marquardt）對垂直漸近線非常敏感，因為 x 中的小誤差會被大大放大.

誰能指出我可以解決這個問題的合適策略？

我有一些統計學的基本知識，但這仍然非常有限。我會渴望學習，只要我知道從哪裡開始尋找:)

非常感謝您的建議！

編輯請原諒忘記提及錯誤。唯一明顯的噪音在，而且是加法的。

編輯 2有關此問題背景的一些附加信息。上圖模擬了聚合物的拉伸行為。正如@whuber 在評論中指出的那樣，您需要得到一個像上面的圖表。

至於到目前為止人們是如何擬合這條曲線的：似乎人們通常會切斷垂直漸近線，直到找到合適的擬合。但是，截止選擇仍然是任意的，這使得擬合過程不可靠且不可重複。

編輯 3&4固定圖。

我們用來手動擬合的方法（即探索性數據分析）可以非常好地處理這些數據。

我希望稍微重新參數化模型以使其參數為正：

對於給定的，讓我們假設有一個獨特的真實滿足這個方程；叫這個或者，為簡潔起見，什麼時候被理解。

我們觀察到一組有序對在哪裡背離通過具有零均值的獨立隨機變量。在本次討論中，我將假設它們都有一個共同的方差，但是這些結果的擴展（使用加權最小二乘法）是可能的、顯而易見的並且易於實現。這是此類集合的模擬示例值，與,, 和一個共同的方差.

**這是一個（故意）艱難的例子，**正如非物質的（消極的）可以理解的那樣值及其非凡的傳播（通常是水平單位，但范圍可達或者在軸）。如果我們能夠對這些數據進行合理擬合，就可以接近估計,，和使用，我們確實會做得很好。

探索性擬合是迭代的。 每個階段包括兩個步驟：估計（基於數據和之前的估計和的和，從中先前的預測值可以獲得) 然後估計. 因為誤差在x中，所以擬合估計來自，而不是反過來。首先在錯誤中排序，什麼時候足夠大，

因此，我們可能會更新通過用最小二乘法擬合這個模型（注意它只有一個參數——斜率，–並且沒有截距）並將係數的倒數作為更新的估計.

接下來，當足夠*小，*反二次項占主導地位，我們發現（再次是誤差中的一階）

再次使用最小二乘法（僅使用斜率項) 我們得到一個更新的估計通過擬合斜率的平方根。

**要了解為什麼這樣有效，**可以通過繪製粗略的探索性近似來獲得這種擬合反對對於較小的. 更好的是，因為測量誤差和單調變化，我們應該關注具有較大值的數據. 這是我們模擬數據集中的一個示例，顯示了最大的一半紅色，最小的一半藍色，一條穿過原點的線適合紅點。

這些點大致對齊，儘管在較小的值處有一點曲率和. （注意軸的選擇：因為是測量值，通常將其繪製在垂直軸上。）通過將擬合集中在曲率應該最小的紅點上，我們應該獲得合理的估計. 的價值標題中顯示的是這條線斜率的平方根：它只是比真實值少%！

此時可以通過以下方式更新預測值

迭代直到估計穩定（不能保證）或者它們在小範圍的值中循環（仍然不能保證）。

事實證明除非我們有一組很好的非常大的值，否則很難估計, 但那- 它確定原始圖中（在問題中）的垂直漸近線並且是問題的焦點 - 可以非常準確地確定，前提是垂直漸近線上有一些數據。 在我們的運行示例中，迭代確實收斂到（這幾乎是正確值的兩倍）和（接近正確的值）。該圖再次顯示了數據，其上疊加了 (a) 灰色的真實曲線（虛線）和 (b) 紅色的估計曲線（實線）：

**這種擬合非常好，以至於很難區分真實曲線和擬合曲線：**它們幾乎在任何地方都重疊。順便說一下，估計的誤差方差為非常接近真實值.

這種方法存在一些問題：

估計是有偏差的。當數據集較小且靠近 x 軸的值相對較少時，偏差變得明顯。適合度系統地有點低。

估計過程需要一種方法來區分“大”和“小”值. 我可以提出探索性的方法來確定最佳定義，但實際上，您可以將這些保留為“調整”常量並更改它們以檢查結果的敏感性。我通過根據值將數據分成三個相等的組來任意設置它們並使用兩個外部組。

該程序不適用於所有可能的組合和或所有可能的數據范圍。但是，只要在數據集中表示了足夠多的曲線以反映兩種漸近線，它就應該能很好地工作：一端是垂直的，另一端是傾斜的。

代碼

以下是用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, xHat} = NestWhile[estimate[##, data] &, {0, 0, data[[1]]}, 
               Norm[Most[#1] - Most[#2]] >= 0.001 &,  2, 100]