python 中使用 curve_fit() 的自定义函数出现 nans 和 infs 错误

Question

我试图找到如何使用自定义函数从python数组中获取高阶矩，例如高斯-埃尔米特函数的h3矩。我将解释我的完整代码设置，以帮助您更好地理解这个问题。

我正在创建星系模拟速度值的统计图（如标准差、h3、h4 等）。我有一个星系模拟快照，其中包含恒星的空间“x_new”、“y_new”值以及占据这些空间坐标的恒星的速度“vz_new”值。我将星系的 xy 平面划分为多个单元，然后让恒星的速度值落入相应的单元中。然后，我通常只使用 np.std()（这是一个内置的 numpy 函数）获取每个单元格中速度值数组的标准差，并使用这些数字绘制地图。

然而，为了获得 h3 或 h4 值，我需要首先从速度值的每个像元中获取频谱，然后定义 Gauss-Hermite 函数，然后获得速度值与该 Gauss-Hermite 函数的曲线拟合，然后获得 h3 或该曲线拟合的 h4 值。因此，我尝试定义自己的函数 get_h3()。这是我遇到一些“NAN”错误的部分，例如：ValueError：数组不得包含 infs 或 NaN 这是我正在使用的代码的相关部分。 ：

import pynbody
import numpy as np
import matplotlib.pyplot as plt
import sys
from scipy import stats
from scipy.stats import norm
import pylab as plb
from scipy.optimize import curve_fit
# import plot lib
import seaborn as sns

#the galaxy spatial coordinates are stored in arrays x_new, y_new. The velocity values are stored in array v_new.


#laying a grid on the galaxy snapshot on an x-y plane
x_new = x_new // 2 + 20
y_new = y_new // 2 + 20
k=(40)
cells = [[[] for y_new in range(k)] for x_new in range(k)] #creating cells or pixels on x-y plane


#letting vz_new values to fall into the grid cells
for ycell in range(k):
    for xcell in range(k):
        cells[ycell][xcell] = vz_new[(y_new  == ycell) & (x_new  == xcell)]
          
for ycell in range(k):
     for xcell in range(k):
        this = cells[ycell][xcell]
      #print(ycell, xcell, len(this), this, sep='\t')
   


 #Define function to get h3 value from velocity arrays in every cell. This will involve nested functions. 
#We will take the velocity array and create a spectrum using plt.hist to get velocity values and corresponding frequencies.

def get_h3(dataV):
    counts, binz, bars = plt.hist(dataV, density=True, bins=150)
    x_h = binz[:-1]
    y_h = counts

    n = len(x_h)                            
    mean = sum(x_h * y_h) / n                   
    sigma = sum(y_h * (x_h - mean)**2) / n        

    
    # Let's define the Gauss-Hermite function
    #a=amplitude, x0=location of peak, sig=std dev, h3, h4
    def gh_func(x_h, a, x0, sig, h3, h4):
        return (
            a * np.exp(-0.5 * ( (x_h - x0) / sig )**2) *
            (
                1 +
                h3 * (
                    ((x_h - x0) / sig)**1 + (-np.sqrt(3)) +
                    ((x_h - x0) / sig)**3 * (2 / np.sqrt(3))
                ) +
                h4 * (
                    (np.sqrt(6) / 4) +
                    ((x_h - x0) / sig)**2 * (-np.sqrt(6)) +
                    ((x_h - x0) / sig)**4 * (np.sqrt(6) / 3)
                )
            )
        )

    popt, pcov = curve_fit(
        gh_func, x_h, y_h, 
        p0=( np.max(y_h), mean, sigma, 0, 0 ), maxfev=5000
    ) #i'm taking h3=0 and h4=0 as guesses
    return popt[3] #this popt[3] is our h3 value

####-----------------------------------------------------------------------------------------
#getting h3 from velocity values in each cell-I THINK THIS IS WHERE THE MAIN ISSUE IS.
h3 = [[[] for y_new in range(k)] for x_new in range(k)]
#cells[ycell][xcell]=
for ycell in range(k):
    for xcell in range(k):
        #cells[ycell][xcell] = np.array(vz_new[(y_new == ycell) & (x_new == xcell)],dtype=object)
        cells[ycell][xcell] = vz_new[(y_new == ycell) & (x_new == xcell)]
        # QUESTION: is there a way to replace invalid nans or infs with 0 for this next part?
        h3[ycell][xcell] = get_h3(cells[ycell][xcell])

fig, ax = plt.subplots()
# heatmap
ax=sns.heatmap(h3, cmap='RdBu')
ax.invert_yaxis()
#plt.savefig('spec_heatmap_150bins_15deg_trial7.png')
plt.show()

有人可以帮助我了解如何处理这些错误或定义一个可以忽略 nan 和 infs 值的函数吗？我的编程经验并不丰富，所以如果有任何不清楚的地方，我深表歉意。请告诉我，我将进一步澄清。多谢。