2016-06-14 11 views
5

散布図の境界にカーブを合わせようとしています。 See this image for reference散布図の境界にカーブを合わせる

私は次の(簡略化された)コードで既に適合しています。これは、データフレームを小さな垂直ストリップにスライスし、幅が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

私はこの問題を本当に興味深いと見出したので、試してみることにしました。私はピジョンまたは自然については知らないが、の情報をすべてポイントで使用しながら、あなたのようなデータセットにエッジをフィッティングするより正確な方法を見つけたと思う。

まず、表示されているようなランダムなデータを生成しましょう。この部分は簡単にスキップすることができます。コードを完全に再現できるように簡単に投稿しています。私は2つの二変量正規分布を使ってそれらの過密をシミュレートし、一様に分布するランダムな点の層を散布しました。その後、彼らはあなたに似たライン式に加え、線の下のすべてが最終的な結果はこのように見ていると、切断されました: 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のような一般的に使用される回帰方法は、ymodel(x)の間の残差が最小になるようにデータ値yを取り、モデルの自由パラメータを最適化します。非線形最小二乗法は、すべてのステップでフィットを改善するために、各ステップでカーブパラメータを揺らぎようとする反復プロセスです。さて、明らかに、これは私たちができるだけ最適なカーブから遠くに私たちを取るように私たちの最小化ルーチンを(あまり遠く)したいので、はしたくないです。

代わりに、次の機能を考慮してください。残差を単純に返すのではなく、反復のすべてのステップで曲線の上のポイントを「反転」し、それらも考慮に入れます。このようにして、カーブの下にあるポイントより効果的に常にポイントが上にあり、繰り返しごとにカーブが下にシフトします。一番低い点に達すると、関数の最小値が見つけられ、散布の端も検出されました。もちろん、この方法では、曲線の下に外れ値がないと仮定しています。しかし、あなたの人物はそれほど苦しんでいないようです。ここで

はこのアイデアを実装する関数です:

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の平均を取るために、より堅牢であることだと思う私も

タックル(およびPythonでそれを実装)しようとしているという、興味深い問題であり、 (または問題に応じてk-最も高い、問題に応じて)データポイントを計算し、平均値を適合させます(適合パラメータが頑強であることを確認する必要があります)k。 このアイデアは、例えばsupplementsPNAS paperにあります。