Tôi thấy vấn đề thực sự thú vị, vì vậy tôi đã quyết định thử. Tôi không biết về pythonic hoặc tự nhiên, nhưng tôi nghĩ rằng tôi đã tìm thấy một cách chính xác hơn để gắn một cạnh vào tập dữ liệu giống như của bạn trong khi sử dụng thông tin từ mọi điểm.
Trước hết, hãy tạo một dữ liệu ngẫu nhiên giống như dữ liệu bạn đã hiển thị. Phần này có thể dễ dàng bỏ qua, tôi đăng nó đơn giản để mã sẽ được hoàn thành và tái sản xuất. Tôi đã sử dụng hai bản phân phối bình thường bivariate để mô phỏng những quá trình đó và rắc chúng với một lớp các điểm ngẫu nhiên được phân bố đều. Sau đó, họ đã được thêm vào một phương trình dòng tương tự như của bạn, và tất cả mọi thứ dưới dòng bị cắt đứt, với kết quả cuối cùng nhìn như thế này: 
Dưới đây là đoạn mã để làm cho nó:
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)
Bây giờ chúng ta có dữ liệu và mô hình, chúng ta có thể suy nghĩ làm thế nào để phù hợp với một cạnh của sự phân bố điểm. Các phương thức hồi quy thường được sử dụng như các ô vuông nhỏ nhất phi tuyến scipy.optimize.curve_fit lấy các giá trị dữ liệu y và tối ưu hóa các tham số miễn phí của một mô hình sao cho số dư giữa y và model(x) là tối thiểu. Các ô vuông nhỏ nhất phi tuyến là một quá trình lặp đi lặp lại cố gắng làm lung lay các tham số đường cong ở mọi bước để cải thiện sự phù hợp ở mọi bước. Bây giờ rõ ràng, đây là một điều chúng tôi không muốn làm, vì chúng tôi muốn quy trình giảm thiểu của chúng tôi đưa chúng tôi đến đường cong phù hợp nhất có thể (nhưng không phải quá ở xa).
Vì vậy, thay vào đó, hãy xem xét chức năng sau. Thay vì chỉ trả lại số dư, nó cũng sẽ "lật" các điểm trên đường cong ở mọi bước của quá trình lặp lại và nhân tố chúng trong đó. Bằng cách này, luôn luôn có nhiều điểm hơn bên dưới đường cong ở trên nó, làm cho đường cong bị dịch chuyển xuống với mỗi lần lặp! Khi đạt đến các điểm thấp nhất, mức tối thiểu của hàm được tìm thấy, và do đó là cạnh của tán xạ. Tất nhiên, phương pháp này giả định bạn không có ngoại lệ bên dưới đường cong - nhưng sau đó con số của bạn dường như không bị nhiều.
Sau đây là các chức năng thực hiện ý tưởng nà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)
Hãy xem cách này có vẻ cho các dữ liệu trên:
# plotting the mock data
plt.plot(dset[0], dset[1], '.', alpha=0.2, label = 'Test data')
# mask bad data (we accidentaly generated some NaN values)
gmask = np.isfinite(dset[1])
dset = dset[np.vstack([gmask, gmask])].reshape((2, -1))
from scipy.optimize import leastsq
guesses =[100, 100, 0]
fit_pars, flag = leastsq(func = flipped_resid, x0 = guesses,
args = (dset[0], dset[1]))
# plot the fit:
y_fit = model(x_data, *fit_pars)
y_guess = model(x_data, *guesses)
plt.plot(x_data, y_fit, 'r-', zorder = 0.9, label = 'Edge')
plt.plot(x_data, y_guess, 'g-', zorder = 0.9, label = 'Guess')
plt.legend(loc = 'lower left')
plt.show()
Phần quan trọng nhất trên là lời kêu gọi leastsq chức năng.Hãy chắc chắn rằng bạn cẩn thận với các dự đoán ban đầu - nếu đoán không rơi vào phân tán, mô hình có thể không hội tụ đúng cách. Sau khi đưa ra dự đoán phù hợp trong ...
Voilà! Các cạnh hoàn toàn phù hợp với thực tế.
Phù hợp với đường cong đến ranh giới của một phân tán