私はこの問題を本当に興味深いと見出したので、試してみることにしました。私はピジョンまたは自然については知らないが、の情報をすべてポイントで使用しながら、あなたのようなデータセットにエッジをフィッティングするより正確な方法を見つけたと思う。
まず、表示されているようなランダムなデータを生成しましょう。この部分は簡単にスキップすることができます。コードを完全に再現できるように簡単に投稿しています。私は2つの二変量正規分布を使ってそれらの過密をシミュレートし、一様に分布するランダムな点の層を散布しました。その後、彼らはあなたに似たライン式に加え、線の下のすべてが最終的な結果はこのように見ていると、切断されました: 
ここでそれを作るためのコードスニペットです:
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
とmodel(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
関数の呼び出しです。あなたが最初の推測に注意していることを確認してください。推測が散布に当たっていなければ、モデルは適切に収束しないかもしれません。で...

ほら、適切な推測を入れた後に!エッジは実際のエッジと完全に一致しています。