如何為非線性最小二乘擬合選擇初始值
上面的問題說明了一切。基本上我的問題是一個通用的擬合函數(可能是任意複雜的),它在我試圖估計的參數中是非線性的,如何選擇初始值來初始化擬合?我正在嘗試做非線性最小二乘。有什麼策略或方法嗎?這個有研究過嗎?有參考嗎?除了臨時猜測之外還有什麼?具體來說,現在我正在使用的擬合形式之一是高斯加線性形式,其中包含我試圖估計的五個參數,例如
在哪裡(橫坐標數據)和(縱坐標數據)意味著在對數空間中我的數據看起來像一條直線加上一個我用高斯近似的凹凸。我沒有理論,沒有什麼可以指導我如何初始化非線性擬合,除了可能像線條的斜率和凹凸的中心/寬度那樣的圖形和眼球觀察。但是我有超過一百個這樣的擬合,而不是圖形和猜測,我更喜歡一些可以自動化的方法。
我在圖書館或網上找不到任何參考資料。我唯一能想到的就是隨機選擇初始值。MATLAB 提供從均勻分佈的 [0,1] 中隨機選擇值。所以對於每個數據集,我運行隨機初始化的擬合一千次,然後選擇最高的? 還有其他(更好的)想法嗎?
附錄#1
首先,這裡有一些數據集的可視化表示,只是為了向你們展示我在說什麼類型的數據。我以原始形式發布數據,沒有任何類型的轉換,然後在日誌空間中發布它的可視化表示,因為它澄清了一些數據的特徵,同時扭曲了其他特徵。我正在發布好數據和壞數據的樣本。
每個圖中的六個面板中的每一個都顯示了四個數據集,紅色、綠色、藍色和青色繪製在一起,每個數據集正好有 20 個數據點。由於數據中看到的顛簸,我試圖用一條直線加上一個高斯來擬合它們中的每一個。
第一個數字是一些好的數據。第二個圖是圖一中相同好的數據的對數圖。第三個數字是一些不良數據。第四個圖是圖三的對數圖。還有更多數據,這些只是兩個子集。大部分數據(大約 3/4)是好的,類似於我在這裡展示的好的數據。
現在有一些評論,請耐心等待,因為這可能會很長,但我認為所有這些細節都是必要的。我會盡量簡潔。
我最初期望一個簡單的冪律(意味著對數空間中的直線)。當我在對數空間中繪製所有內容時,我看到了 4.8 mHz 左右的意外顛簸。顛簸經過徹底調查,並且在其他工作中也被發現,所以不是我們搞砸了。它在物理上就在那裡,其他已發表的作品也提到了這一點。所以我只是在我的線性形式中添加了一個高斯項。請注意,這種擬合是在日誌空間中完成的(因此我的兩個問題包括這個)。
現在,在閱讀了Stumpy Joe Pete 對我的另一個問題的回答(根本與這些數據無關)並閱讀了這個和這個以及其中的引用(Clauset 的東西)之後,我意識到我不應該適合 log-log空間。所以現在我想在預先轉換的空間中做所有事情。
問題一:看好數據,我還是覺得在預變換空間中一個線性加一個高斯還是一個好形式。我很想听聽其他擁有更多數據經驗的人的想法。高斯+線性合理嗎?我應該只做高斯嗎?還是完全不同的形式?
問題 2:無論問題 1 的答案如何,我仍然需要(很可能)非線性最小二乘擬合,因此仍然需要初始化幫助。
我們看到兩組數據,我們非常喜歡在 4-5 mHz 左右捕獲第一個顛簸。所以我不想添加更多的高斯項,我們的高斯項應該以第一個凹凸為中心,這幾乎總是更大的凹凸。我們希望在 0.8mHz 和 5mHz 左右之間“更準確”。我們不太關心更高的頻率,但也不想完全忽略它們。所以也許某種稱重?或者 B 可以始終在 4.8mHz 左右初始化?
橫坐標數據是以毫赫茲為單位的頻率,表示為. 縱坐標數據是我們正在計算的係數,表示為. 所以沒有對數轉換,形式為
- 是頻率,總是積極的。
- 是正係數。所以我們在第一象限工作。
- ,我認為幅度也應該始終為正,因為我們只是在處理顛簸。當我查看數據時,我總是看到高峰而沒有低谷。看起來在所有數據中都有多個更高頻率的顛簸。第一個凸起總是比其他凸起大得多。在良好的數據中,次級凸起非常弱,但在不良數據(例如面板 2 和 5)中,次級凸起很強。所以我們實際上沒有一個山谷,而是兩個顛簸。意味著幅度. 而且由於我們主要關心第一個高峰,因此更有理由要積極。
- 是凹凸的中心,我們總是希望它在 4-5mHz 左右的那個大凹凸上。在我們解析的頻率範圍內,它幾乎總是出現在 4.8mHz。
- 是凹凸的寬度。我想它是關於零意義對稱的將具有相同的效果因為它是平方的。所以我們不在乎它的價值是什麼。假設我們更喜歡它是積極的。
- 是線的斜率,似乎它可能略微為負,因此不對它實施任何限制。斜坡本身就很有趣,所以我們不想對其施加任何限制,我們只想看看它會是什麼。它是積極的還是消極的?它的大小有多大/小?等等。
- 是(幾乎)-截距。這裡的微妙之處在於,由於高斯項,不完全是-截距。實際截距(如果我們推斷) 將是
所以這裡唯一的限制是截距也應該是正數。截距為零,我不知道這意味著什麼。但消極肯定看起來很荒謬。我想在這裡我們可以允許如有必要,可略為負,幅度較小。原因截距在這裡很重要,但我們的一些同事實際上也對外推感興趣。我們擁有的最小頻率是 0.8mHz,他們想在 0 到 0.8mHz 之間進行推斷。我幼稚的想法是只使用適合的方式一直下降到.
我知道外推比插值更難/更危險,但使用直線加高斯(希望它衰減得足夠快)對我來說似乎是合理的。有點像具有自然邊界條件的自然三次樣條曲線,左端點處的斜率,只需延長線並查看它與軸。如果它不是負數,則使用該線進行外推。
問題3:你們認為在這種情況下以這種方式推斷是什麼?有什麼優點/缺點嗎?還有其他外推想法嗎?同樣,我們只關心較低的頻率,因此在 0 和 1mHz 之間進行外推……有時非常非常小的頻率,接近於零。我知道這篇文章已經打包了。我在這裡問了這個問題,因為答案可能是相關的,但如果你們願意,我可以把這個問題分開,稍後再問另一個。
最後,根據要求,這裡有兩個示例數據集。
0.813010000000000 0.091178000000000 0.012728000000000 1.626000000000000 0.103120000000000 0.019204000000000 2.439000000000000 0.114060000000000 0.063494000000000 3.252000000000000 0.123130000000000 0.071107000000000 4.065000000000000 0.128540000000000 0.073293000000000 4.878000000000000 0.137040000000000 0.074329000000000 5.691100000000000 0.124660000000000 0.071992000000000 6.504099999999999 0.104480000000000 0.071463000000000 7.317100000000000 0.088040000000000 0.070336000000000 8.130099999999999 0.080532000000000 0.036453000000000 8.943100000000001 0.070902000000000 0.024649000000000 9.756100000000000 0.061444000000000 0.024397000000000 10.569000000000001 0.056583000000000 0.025222000000000 11.382000000000000 0.052836000000000 0.024576000000000 12.194999999999999 0.048727000000000 0.026598000000000 13.008000000000001 0.045870000000000 0.029321000000000 13.821000000000000 0.041454000000000 0.067300000000000 14.633999999999999 0.039596000000000 0.081800000000000 15.447000000000001 0.038365000000000 0.076443000000000 16.260000000000002 0.036425000000000 0.075912000000000
第一列是以 mHz 為單位的頻率,在每個數據集中都相同。第二列是一個好的數據集(好的數據圖一和二,面板 5,紅色標記),第三列是一個壞的數據集(壞數據圖三和四,面板 5,紅色標記)。
希望這足以激發一些更開明的討論。謝謝大家。
如果有一種既好又通用的策略——總是有效的——它已經在每個非線性最小二乘程序中實施,並且起始值將不是問題。
對於許多特定問題或問題系列,存在一些很好的初始值方法;一些軟件包帶有針對特定非線性模型的良好起始值計算或更通用的方法,這些方法通常有效,但可能需要通過更具體的函數或直接輸入起始值來幫助解決。
在某些情況下,探索這個空間是必要的,但我認為你的情況可能是這樣的,更具體的策略可能是值得的——但設計一個好的策略幾乎需要很多我們不太可能擁有的領域知識。
對於您的特定問題,可以設計可能的好策略,但這不是一個簡單的過程;您對峰的可能大小和範圍(典型參數值和典型’s會給出一些想法),可以做更多的事情來設計一個好的起始值選擇。
什麼是典型的範圍’沙你明白了嗎?平均結果是什麼樣的?有哪些不尋常的案例?您知道哪些參數值是可能的或不可能的?
一個例子——高斯部分是否一定會產生一個轉折點(峰或谷)?或者它有時相對於線太小以至於它沒有?是總是積極的?等等。
如果可以的話,一些示例數據會有所幫助 - 典型案例和困難案例。
編輯:這是一個例子,如果問題不是太嘈雜,你可以如何做得相當好:
這是從您的模型生成的一些數據(總體值為 A = 1.9947,B = 10,C = 2.828,D = 0.09,E = 5):
我能夠估計的起始值是
(As = 1.658,Bs = 10.001,Cs = 3.053,Ds = 0.0881,Es = 5.026)
該啟動模型的擬合如下所示:
步驟是:
- 擬合 Theil 回歸以獲得 D 和 E 的粗略估計
- 減去 Theil 回歸的擬合
- 使用 LOESS 擬合平滑曲線
- 找到峰值得到A的粗略估計,峰值對應的x值得到B的粗略估計
- 將 y 值大於 A 估計值的 60% 的 LOESS 擬合作為觀測值並擬合二次
- 使用二次來更新 B 的估計並估計 C
- 從原始數據中減去高斯的估計值
- 將另一個 Theil 回歸擬合到該調整後的數據以更新 D 和 E 的估計值
在這種情況下,這些值非常適合開始非線性擬合。
我把它寫成
R
代碼,但同樣的事情可以在 MATLAB 中完成。我認為比這更好的事情是可能的。
如果數據非常嘈雜,這根本不會奏效。
Edit2:這是我在 R 中使用的代碼,如果有人感興趣的話:
gausslin.start <- function(x,y) { theilreg <- function(x,y){ yy <- outer(y, y, "-") xx <- outer(x, x, "-") z <- yy / xx slope <- median(z[lower.tri(z)]) intercept <- median(y - slope * x) cbind(intercept=intercept,slope=slope) } tr <- theilreg(x,y1) abline(tr,col=4) Ds = tr[2] Es = tr[1] yf <- y1-Ds*x-Es yfl <- loess(yf~x,span=.5) # assumes there are enough points that the maximum there is 'close enough' to # the true maximum yflf <- yfl$fitted locmax <- yflf==max(yflf) Bs <- x[locmax] As <- yflf[locmax] qs <- yflf>.6*As ys <- yfl$fitted[qs] xs <- x[qs]-Bs lf <- lm(ys~xs+I(xs^2)) bets <- lf$coefficients Bso <- Bs Bs <- Bso-bets[2]/bets[3]/2 Cs <- sqrt(-1/bets[3]) ystart <- As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es y1a <- y1-As*exp(-((x-Bs)/Cs)^2) tr <- theilreg(x,y1a) Ds <- tr[2] Es <- tr[1] res <- data.frame(As=As, Bs=Bs, Cs=Cs, Ds=Ds, Es=Es) res }
.
# population parameters: A = 1.9947 , B = 10, C = 2.828, D = 0.09, E = 5 # generate some data set.seed(seed=3424921) x <- runif(50,1,30) y <- dnorm(x,10,2)*10+rnorm(50,0,.2) y1 <- y+5+x*.09 # This is the data xo <- order(x) starts <- gausslin.start(x,y1) ystart <- with(starts, As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es) plot(x,y1) lines(x[xo],ystart[xo],col=2)