python 三维曲线拟合_基于三维数据和参数的Scipy曲线拟合-CSDN博客

我正致力于在scipy中拟合三维分布函数。我有一个numpy数组，在x和y-bin中有计数，我正试图将其与一个相当复杂的三维分布函数相匹配。数据适合26(！)描述其两个组成种群形状的参数。

我在这里了解到，当我调用leatsq时，必须将x和y坐标作为“args”传递。unutbu提供的代码是为我编写的，但是当我试图将其应用于我的特定情况时，会出现错误“TypeError:leastsq()为关键字参数“args”获取多个值”

这是我的代码(对不起，长度太长了)：import numpy as np

import matplotlib.pyplot as plt

import scipy.optimize as spopt

from textwrap import wrap

import collections

cl = 0.5

ch = 3.5

rl = -23.5

rh = -18.5

mbins = 10

cbins = 10

def hist_data(mixed_data, mbins, cbins):

import numpy as np

H, xedges, yedges = np.histogram2d(mixed_data[:,1], mixed_data[:,2], bins = (mbins, cbins), weights = mixed_data[:,3])

x, y = 0.5 * (xedges[:-1] + xedges[1:]), 0.5 * (yedges[:-1] + yedges[1:])

return H.T, x, y

def gauss(x, s, mu, a):

import numpy as np

return a * np.exp(-((x - mu)**2. / (2. * s**2.)))

def tanhlin(x, p0, p1, q0, q1, q2):

import numpy as np

return p0 + p1 * (x + 20.) + q0 * np.tanh((x - q1)/q2)

def func3d(p, x, y):

import numpy as np

from sys import exit

rsp0, rsp1, rsq0, rsq1, rsq2, rmp0, rmp1, rmq0, rmq1, rmq2, rs, rm, ra, bsp0, bsp1, bsq0, bsq1, bsq2, bmp0, bmp1, bmq0, bmq1, bmq2, bs, bm, ba = p

x, y = np.meshgrid(coords[0], coords[1])

rs = tanhlin(x, rsp0, rsp1, rsq0, rsq1, rsq2)

rm = tanhlin(x, rmp0, rmp1, rmq0, rmq1, rmq2)

ra = schechter(x, rap, raa, ram) # unused

bs = tanhlin(x, bsp0, bsp1, bsq0, bsq1, bsq2)

bm = tanhlin(x, bmp0, bmp1, bmq0, bmq1, bmq2)

ba = schechter(x, bap, baa, bam) # unused

red_dist = ra / (rs * np.sqrt(2 * np.pi)) * gauss(y, rs, rm, ra)

blue_dist = ba / (bs * np.sqrt(2 * np.pi)) * gauss(y, bs, bm, ba)

result = red_dist + blue_dist

return result

def residual(p, coords, data):

import numpy as np

model = func3d(p, coords)

res = (model.flatten() - data.flatten())

# can put parameter restrictions in here

return res

def poiss_err(data):

import numpy as np

return np.where(np.sqrt(H) > 0., np.sqrt(H), 2.)

# =====

H, x, y = hist_data(mixed_data, mbins, cbins)

data = H

coords = x, y

# x and y will be the projected coordinates of the data H onto the plane z = 0

# x has bins of width 0.5, with centers at -23.25, -22.75, ... , -19.25, -18.75

# y has bins of width 0.3, with centers at 0.65, 0.95, ... , 3.05, 3.35

Param = collections.namedtuple('Param', 'rsp0 rsp1 rsq0 rsq1 rsq2 rmp0 rmp1 rmq0 rmq1 rmq2 rs rm ra bsp0 bsp1 bsq0 bsq1 bsq2 bmp0 bmp1 bmq0 bmq1 bmq2 bs bm ba')

p_guess = Param(rsp0 = 0.152, rsp1 = 0.008, rsq0 = 0.044, rsq1 = -19.91, rsq2 = 0.94, rmp0 = 2.279, rmp1 = -0.037, rmq0 = -0.108, rmq1 = -19.81, rmq2 = 0.96, rs = 1., rm = -20.5, ra = 10000., bsp0 = 0.298, bsp1 = 0.014, bsq0 = -0.067, bsq1 = -19.90, bsq2 = 0.58, bmp0 = 1.790, bmp1 = -0.053, bmq0 = -0.363, bmq1 = -20.75, bmq2 = 1.12, bs = 1., bm = -20., ba = 2000.)

opt, cov, infodict, mesg, ier = spopt.leastsq(residual, p_guess, poiss_err(H), args = coords, maxfev = 100000, full_output = True)

这是我的数据，只有更少的箱子：[[ 1.00000000e+01 1.10000000e+01 2.10000000e+01 1.90000000e+01

1.70000000e+01 2.10000000e+01 2.40000000e+01 1.90000000e+01

2.80000000e+01 1.90000000e+01]

[ 1.40000000e+01 4.50000000e+01 6.00000000e+01 6.80000000e+01

1.34000000e+02 1.97000000e+02 2.23000000e+02 2.90000000e+02

3.23000000e+02 3.03000000e+02]

[ 3.00000000e+01 1.17000000e+02 3.78000000e+02 9.74000000e+02

1.71900000e+03 2.27700000e+03 2.39000000e+03 2.25500000e+03

1.85600000e+03 1.31000000e+03]

[ 1.52000000e+02 9.32000000e+02 2.89000000e+03 5.23800000e+03

6.66200000e+03 6.19100000e+03 4.54900000e+03 3.14600000e+03