2016-06-14 129 views
5

我試圖用曲線的散點圖的邊界的邊界。 See this image for referenceenter image description here曲線擬合,以散點圖

我已經完成了以下(簡化)代碼的適配。它將數據幀分成小垂直條,然後在寬度爲width的那些條中找到最小值,而忽略nan。 (該功能是單調遞減的。)

def func(val): 
    """ returns some function of 'val'""" 
    return val * 2 

for i in range(0, max_val, width)): 
    _df = df[(df.val > i) & (df.val < i + width)] # vertical slice 
    if np.isnan(np.min(func(_df.val)):   # ignore nans 
     continue 
    xs.append(i + width)       
    ys.append(np.min(func(_df.val))) 

我那麼做配合scipy.optimize.curve_fit。我的問題是:是否有更自然或pythonic的方式來做到這一點 - 是否有任何方法可以提高準確性? (例如,通過賦予更高的權重與點的高密度散點圖的地方?)

回答

3

我發現這個問題非常有趣,所以我決定試一試。我不知道Pythonic或自然,但我認爲我已經找到了一種更準確的方式來適合像你這樣的數據集邊緣,同時使用點的信息。

首先,讓我們產生看起來像你所示的隨機數據。這部分可以很容易地跳過,我簡單地發佈它,以便代碼完整和可重複。我使用了兩個二元正態分佈來模擬這些過度密度,並將它們分佈在一層均勻分佈的隨機點上。然後,他們被添加到類似你這樣的直線方程,一切線下被切斷,與最終的結果看起來像這樣: making of the data

下面的代碼片段,使其:

import numpy as np 

x_res = 1000 
x_data = np.linspace(0, 2000, x_res) 

# true parameters and a function that takes them 
true_pars = [80, 70, -5] 
model = lambda x, a, b, c: (a/np.sqrt(x + b) + c) 
y_truth = model(x_data, *true_pars) 

mu_prim, mu_sec = [1750, 0], [450, 1.5] 
cov_prim = [[300**2, 0  ], 
      [  0, 0.2**2]] 
# covariance matrix of the second dist is trickier 
cov_sec = [[200**2, -1  ], 
      [ -1, 1.0**2]] 
prim = np.random.multivariate_normal(mu_prim, cov_prim, x_res*10).T 
sec = np.random.multivariate_normal(mu_sec, cov_sec, x_res*1).T 
uni = np.vstack([x_data, np.random.rand(x_res) * 7]) 

# censoring points that will end up below the curve 
prim = prim[np.vstack([[prim[1] > 0], [prim[1] > 0]])].reshape(2, -1) 
sec = sec[np.vstack([[sec[1] > 0], [sec[1] > 0]])].reshape(2, -1) 

# rescaling to data 
for dset in [uni, sec, prim]: 
    dset[1] += model(dset[0], *true_pars) 

# this code block generates the figure above: 
import matplotlib.pylab as plt 
plt.figure() 
plt.plot(prim[0], prim[1], '.', alpha=0.1, label = '2D Gaussian #1') 
plt.plot(sec[0], sec[1], '.', alpha=0.5, label = '2D Gaussian #2') 
plt.plot(uni[0], uni[1], '.', alpha=0.5, label = 'Uniform') 
plt.plot(x_data, y_truth, 'k:', lw = 3, zorder = 1.0, label = 'True edge') 
plt.xlim(0, 2000) 
plt.ylim(-8, 6) 
plt.legend(loc = 'lower left') 
plt.show() 

# mashing it all together 
dset = np.concatenate([prim, sec, uni], axis = 1) 

現在我們已經有了數據和模型,我們可以集體討論如何擬合點分佈的邊緣。常用的迴歸方法如非線性最小二乘scipy.optimize.curve_fit取數據值y和優化模型的自由參數,使得ymodel(x)之間的殘餘是最小的。非線性最小二乘是一個迭代過程,它試圖在每一步中擺動曲線參數以改善每一步的擬合。現在清楚了,這是一兩件事,我們想做的事,因爲我們希望最小化過程拿我們當遠離最佳擬合曲線越好(但不能太遠)。

所以取而代之,讓我們考慮下面的函數。它不是簡單地返回殘差,而是在迭代的每一步都「翻轉」曲線上方的點,並將它們也考慮在內。通過這種方式,曲線下方的曲線總是比曲線下方的曲線更多,從而導致曲線每次迭代都向下移動!一旦達到最低點,就會發現函數的最小值,散射的邊緣也是如此。當然,這種方法假定你沒有低於曲線的異常值,但是你的數字似乎並沒有受到太多的影響。

下面是實現這個想法的功能:

def get_flipped(y_data, y_model): 
    flipped = y_model - y_data 
    flipped[flipped > 0] = 0 
    return flipped 

def flipped_resid(pars, x, y): 
    """ 
    For every iteration, everything above the currently proposed 
    curve is going to be mirrored down, so that the next iterations 
    is going to progressively shift downwards. 
    """ 
    y_model = model(x, *pars) 
    flipped = get_flipped(y, y_model) 
    resid = np.square(y + flipped - y_model) 
    #print pars, resid.sum() # uncomment to check the iteration parameters 
    return np.nan_to_num(resid) 

讓我們來看看如何查找上面的數據:

​​3210

上面最重要的部分是調用leastsq功能。確保你對最初的猜測小心 - 如果猜測不落在散射點上,那麼模型可能不會收斂。把一個適當的猜測中......

ifitfitsisits

瞧了!邊緣與真實的邊緣完美匹配。

0

這是一個有趣的問題,那我也試圖解決(在蟒蛇實現它)

我認爲與其服用min,更健壯採取的k -lowest平均(或k-最高,取決於問題)數據點,並擬合平均值(還應該檢查擬合參數是否穩健)。 您可以在此PNAS papersupplements 中找到該想法。