2d fit из гауссиана не подойдет - PullRequest
Новая
       71

2d fit из гауссиана не подойдет

0 голосов
/ 04 мая 2018

Я пытаюсь установить 2d набор данных. Я могу построить его без каких-либо проблем, но когда я пытаюсь подогнать его, я всегда получаю одно и то же сообщение об ошибке, которое не понимаю: 

def twod_Gaussian(x,y,Amp,x0,y0,sigma_x,sigma_y,Offset):
    return Amp*np.exp(-((x-x0)/(2*sigma_x**2)+(y-y0)/(2*sigma_y**2)))+Offset

# get data
filename = '20180503-1455-43_confocal_xy_data.dat'
data = np.genfromtxt(filename, comments='#')
values = open(filename)
val = values.readlines(500)
res = float(val[11][-4:-1])
values.close()
x = [float(i) for i in data[:, 0]]
y = [float(i) for i in data[:, 1]]
count_rate = [float(i) for i in data[:, 3]]



#reshape data
xmax = max(x)
xmin = min(x)
ymax = max(y)
ymin = min(y)
cmax = max(count_rate)
cmin = min(count_rate)
np.percentile(count_rate, 99)
resX = (xmax-xmin)/res
rexY = (ymax-ymin)/res
res2 = len(x)/res
xx = np.array(x).reshape((int(res2), int(res)))
yy = np.array(y).reshape((int(res2), int(res)))



#plot data
count_matrix = np.array(count_rate).reshape((int(res2), int(res)))
fig = plt.figure()
ax = fig.gca(projection='3d')
CS = ax.plot_surface((xx)/1e-6, (yy)/1e-6, count_matrix, cmap='plasma')
plt.ticklabel_format(style='plain', axis='both', scilimits=(0, 0))
ax.set_xlabel(r'X position ($\mathrm{\mu}$m)')
ax.set_ylabel(r'Y position ($\mathrm{\mu}$m)')
fig.colorbar(CS, ax=ax, extend='max', format='%.0e')
plt.show()

#initial guess and trying to fit it
p0=(98e-06,81.5e-06,50000,0.5e-06,0.5e-06,20000)
popt, pcov = curve_fit(twod_Gaussian, x, y, count_rate, p0)

Я всегда получаю сообщение об ошибке

ValueError: sigma имеет неправильную форму.

Если я изменю последнюю строку на 

popt, pcov = curve_fit(twod_Gaussian, (x, y), count_rate, p0)

Я получаю сообщение об ошибке: 

twod_Gaussian() missing 1 required positional argument: 'Offset'

Понятия не имею, почему форма должна быть неправильной?

1 Ответ

0 голосов
/ 13 мая 2018

Во-первых, я считаю, что проблема в том, что у определения функции модели есть дополнительный аргумент, потому что вы используете два отдельных аргумента для двух измерений (x и y) вместо одного аргумента для независимой переменной, который curve_fit требуется. Изменение функции модели на 

def twod_Gaussian(xy, Amp,x0,y0,sigma_x,sigma_y,Offset):
    x, y = xy
    return Amp*np.exp(-((x-x0)/(2*sigma_x**2)+(y-y0)/(2*sigma_y**2)))+Offset

, а затем позвонить curve_fit с 

xy = x, y
popt, pcov = curve_fit(twod_Gaussian, xy, count_rate, p0)

должен делать то, что вы хотите.

Во-вторых, вы также можете рассмотреть возможность использования lmfit (https://lmfit.github.io/lmfit-py).). Эта библиотека имеет немного более высокий уровень интерфейса для подбора кривой, все еще обернутый вокруг решателей scipy.optimize. Например, она обрабатывает параметры как первые -класс, именованные объекты, и легко поддерживает размещение границ и ограничений для параметров вне функции модели. В вашем случае он также поддерживает определение нескольких независимых переменных для функции подгонки, не заставляя вас предполагать, что это только первый аргумент модели Функция - это независимая переменная, а все остальные должны быть переменными в фитинге. С lmfit ваша проблема выглядела бы примерно так (я не запускал это, потому что у меня нет ваших данных): 

from lmfit import Model
def twod_Gaussian(x, y, Amp,x0,y0,sigma_x,sigma_y,Offset):
    return Amp*np.exp(-((x-x0)/(2*sigma_x**2)+(y-y0)/(2*sigma_y**2)))+Offset

# turn your model function into a Model, specifying independent variables
g2model = Model(twod_Gaussian, independent_vars=['x', 'y'])

# create Parameters object -- ordered dict of *named parameters, 
# giving initial values here:
params = g2model.make_params(Amp=98e-6, x0=81e-6, y0=50000., sigma_x=5e-7, sigma_y=5e-7, Offset=20000)

# Note that using names is *way* better than using an ordered list. 
# for example, you may have mixed up the order and really meant:
params = g2model.make_params(x0=98e-6, y0=81e-6, Amp=50000., sigma_x=5e-7, sigma_y=5e-7, Offset=20000)

# you can place bounds on parameters, perhaps as
params['sigma_x'].min = 0
params['sigma_y'].min = 0

# do fit, passing independent vars explicitly
result = g2model.fit(count_rate, params, x=x, y=y)

# print out full report of statistics, best-fit values and stderrs:
print(result.fit_report())

# array of best-fit data == result.best_fit

Есть много других функций, описанных в документации.

Наконец, чтобы сделать двумерный гауссиан с вашей функцией, y должен быть массивом по строкам (x по столбцам). В противном случае ваш twod_Gaussian() вернет 1-й массив той же длины, что и x и y. Итак, я думаю, что вы хотите добавить 

y = np.vstack(y)

сразу после создания y из массива data.

...