Python二维高斯拟合

2024-03-30 00:52:53 发布

您现在位置:Python中文网/ 问答频道 /正文

我对Python很陌生,但我正在尝试为一些数据生成2D高斯拟合。具体地说,恒星通量与坐标系/网格中的某些位置有关。然而,在我的网格中并不是所有的位置都有相应的通量值。我真的不想把这些值设置为零,以防它影响我的拟合度,但我似乎不能将它们设置为nan并且仍然使高斯拟合工作。这是我使用的代码(从here稍作修改):

import numpy
import scipy
from numpy import *
from scipy import optimize

def gaussian(height, center_x, center_y, width_x, width_y):
    width_x = float(width_x)
    width_y = float(width_y)
    return lambda x,y: height*exp(-(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)

def moments(data):
    total = nansum(data)
    X, Y = indices(data.shape)
    center_x = nansum(X*data)/total
    center_y = nansum(Y*data)/total
    row = data[int(center_x), :]
    col = data[:, int(center_y)]
    width_x = nansum(sqrt(abs((arange(col.size)-center_y)**2*col))/nansum(col))
    width_y = nansum(sqrt(abs((arange(row.size)-center_x)**2*row))/nansum(row))
    height = nanmax(data)
    return height, center_x, center_y, width_x, width_y

def fitgaussian(data):
    params = moments(data)
    errorfunction = lambda p: ravel(gaussian(*p)(*indices(data.shape)) - data)
    p, success = optimize.leastsq(errorfunction, params)
    return p

parameters = fitgaussian(data)
fit = gaussian(*parameters)

我的通量值在一个名为data的2D数组中。如果我在这个数组中有0而不是nan值,那么这段代码就可以工作了,否则我的parameters总是以[nan nan nan nan nan]的形式出现。如果有办法解决这个问题,我将非常感谢你的洞察力!解释得越详细越好。提前谢谢!在


Tags: importdatareturndefcolgaussiannanwidth
2条回答

显然要做的是从data中删除nan。然而,这样做还需要移除2D XY位置数组中的相应位置:

X, Y = np.indices(data.shape)
mask = ~np.isnan(data)
x = X[mask]
y = Y[mask]
data = data[mask]

现在您可以使用optimize.leastsq(或更新、更简单的optimize.curve_fit) 要使数据适合模型函数:

^{pr2}$

例如,如果我们用nan生成一些随机的data

data = make_data(shape)

所以

import matplotlib.pyplot as plt
plt.imshow(data)
plt.show()

看起来像

enter image description here

用白点表示有NaN值的地方,然后

import numpy as np
from scipy import optimize
np.set_printoptions(precision=4)


def gaussian(p, x, y):
    height, center_x, center_y, width_x, width_y = p
    return height*np.exp(-(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)

def moments(data):
    total = np.nansum(data)
    X, Y = np.indices(data.shape)
    center_x = np.nansum(X*data)/total
    center_y = np.nansum(Y*data)/total
    row = data[int(center_x), :]
    col = data[:, int(center_y)]
    width_x = np.nansum(np.sqrt(abs((np.arange(col.size)-center_y)**2*col))
                        /np.nansum(col))
    width_y = np.nansum(np.sqrt(abs((np.arange(row.size)-center_x)**2*row))
                        /np.nansum(row))
    height = np.nanmax(data)
    return height, center_x, center_y, width_x, width_y

def errorfunction(p, x, y, data):
    return gaussian(p, x, y) - data

def fitgaussian(data):
    params = moments(data)
    X, Y = np.indices(data.shape)
    mask = ~np.isnan(data)
    x = X[mask]
    y = Y[mask]
    data = data[mask]
    p, success = optimize.leastsq(errorfunction, params, args=(x, y, data))
    return p

def make_data(shape):
    h, w = shape
    p = 50, h/2.0, w/2.0, h/3.0, w/5.0
    print('Actual parameters: {}'.format(np.array(p)))
    X, Y = np.indices(shape)
    data = gaussian(p, X, Y) + np.random.random(shape)
    mask = np.random.random(shape) < 0.3
    data[mask] = np.nan
    return data

shape = 100, 200
data = make_data(shape)
X, Y = np.indices(shape)
parameters = fitgaussian(data)
print('Fitted parameters: {}'.format(parameters))
fit = gaussian(parameters, X, Y)

收益率

Actual parameters: [  50.       50.      100.       33.3333   40.    ]
Fitted parameters: [ 50.2908  49.9992  99.9927  33.7039  40.6149]

只需移除所有没有对应通量值的值。如果在这一点上y轴上没有任何内容,则删除成对的值将无关紧要。在

如果空值等于'',则应该删除所有没有通量值的值

# assumes data.shape = (1, 3) where data[:,0:1] is the x,y axis
# data[:,2] contains the flux values
data = numpy.delete(data, numpy.where(data[:,3] == ''), axis=0)

如果空值等于nan,则此操作将完成

^{pr2}$

相关问题 更多 >