英文:
Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?
问题
协方差未在SciPy的Curvefit中估算
这是我的数据集:
frequency (Hz) brightness (ergs/s/cm^2/sr/Hz) brightness (J/s/m^2/sr/Hz)
float64 float64 float64
34473577711.372055 7.029471536390586e-16 7.029471536390586e-19
42896956937.69582 1.0253178228238486e-15 1.0253178228238486e-18
51322332225.44733 1.3544045476166584e-15 1.3544045476166584e-18
60344529880.18272 1.6902073280174815e-15 1.6902073280174815e-18
68767909106.5062 2.0125779972022745e-15 2.0125779972022743e-18
77780126454.10146 2.3148004995630144e-15 2.3148004995630145e-18
... ... ...
489996752265.52826 3.201319839821188e-16 3.201319839821188e-19
506039097962.6759 2.5968748350997043e-16 2.596874835099704e-19
523273092332.3638 2.0595903864583913e-16 2.0595903864583912e-19
539918248580.7806 1.7237876060575648e-16 1.7237876060575649e-19
557158231134.7507 1.3879848256567381e-16 1.3879848256567383e-19
573803387383.1646 1.0521820452559118e-16 1.0521820452559118e-19
591049358121.42 9.178609330955852e-17 9.178609330955852e-20
我尝试使用CurveFit将其拟合到普朗克的辐射定律:
import numpy as np
from scipy.optimize import curve_fit
h = 6.626 * 10e-34
c = 3 * 10e8
k = 1.38 * 10e-23
const1 = 2 * h / (c**2)
const2 = h / k
def planck(x, v):
return const1 * (v**3) * (1 / ((np.exp(const2 * v / x)) - 1))
popt, pcov = curve_fit(planck, cmb['frequency (Hz)'], cmb['brightness (J/s/m^2/sr/Hz)'])
print(popt, pcov)
警告:
/tmp/ipykernel_2500/4072287013.py:11: RuntimeWarning: divide by zero encountered in divide
return const1*(v**3)*(1/((np.exp((const2)*v/x))-1))
我得到了 popt=1
和 pcov=nan
。现在函数中的指数项相差几个数量级。并且一些值不允许数学上近似这个定律。我尝试使用该定律的对数形式,但那也不起作用。我该如何克服这个问题?
英文:
Covariance not estimated in SciPy's Curvefit
Here's my dataset:
frequency (Hz) brightness (ergs/s/cm^2/sr/Hz) brightness (J/s/m^2/sr/Hz)
float64 float64 float64
34473577711.372055 7.029471536390586e-16 7.029471536390586e-19
42896956937.69582 1.0253178228238486e-15 1.0253178228238486e-18
51322332225.44733 1.3544045476166584e-15 1.3544045476166584e-18
60344529880.18272 1.6902073280174815e-15 1.6902073280174815e-18
68767909106.5062 2.0125779972022745e-15 2.0125779972022743e-18
77780126454.10146 2.3148004995630144e-15 2.3148004995630145e-18
... ... ...
489996752265.52826 3.201319839821188e-16 3.201319839821188e-19
506039097962.6759 2.5968748350997043e-16 2.596874835099704e-19
523273092332.3638 2.0595903864583913e-16 2.0595903864583912e-19
539918248580.7806 1.7237876060575648e-16 1.7237876060575649e-19
557158231134.7507 1.3879848256567381e-16 1.3879848256567383e-19
573803387383.1646 1.0521820452559118e-16 1.0521820452559118e-19
591049358121.42 9.178609330955852e-17 9.178609330955852e-20
I tried to use CurveFit to fit this to Planck's Radiation Law:
import numpy as np
from scipy.optimize import curve_fit
h=6.626*10e-34
c=3*10e8
k=1.38*10e-23
const1=2*h/(c**2)
const2=h/k
def planck(x,v):
return const1*(v**3)*(1/((np.exp(const2*v/x))-1))
popt,pcov= curve_fit(planck, cmb['frequency (Hz)'],cmb['brightness (J/s/m^2/sr/Hz)'])
print(popt, pcov)
Warning:
/tmp/ipykernel_2500/4072287013.py:11: RuntimeWarning: divide by zero encountered in divide
return const1*(v**3)*(1/((np.exp((const2)*v/x))-1))
I get popt=1
and pcov=nan
. Now the exponential term in the function differs by several orders of magnitude. And some of the values don't permit to approximate the law mathematically. I tried using the logarithmic form of the law but that doesn't work either. How can I overcome this problem?
答案1
得分: 1
有很多问题,包括你的变量被交换,你不必要地重新定义物理常数,而且你的表达式在数值上不稳定。你需要使用exp1m
代替:
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import h, c, k
from scipy.optimize import curve_fit
freq, brightness_erg, brightness_j = np.array(((
(34473577711.372055, 7.0294715363905860e-16, 7.0294715363905860e-19),
(42896956937.695820, 1.0253178228238486e-15, 1.0253178228238486e-18),
(51322332225.447330, 1.3544045476166584e-15, 1.3544045476166584e-18),
(60344529880.182720, 1.6902073280174815e-15, 1.6902073280174815e-18),
(68767909106.506200, 2.0125779972022745e-15, 2.0125779972022743e-18),
(77780126454.101460, 2.3148004995630144e-15, 2.3148004995630145e-18),
(489996752265.52826, 3.2013198398211880e-16, 3.2013198398211880e-19),
(506039097962.67590, 2.5968748350997043e-16, 2.5968748350997040e-19),
(523273092332.36380, 2.0595903864583913e-16, 2.0595903864583912e-19),
(539918248580.78060, 1.7237876060575648e-16, 1.7237876060575649e-19),
(557158231134.75070, 1.3879848256567381e-16, 1.3879848256567383e-19),
(573803387383.16460, 1.0521820452559118e-16, 1.0521820452559118e-19),
(591049358121.42000, 9.1786093309558520e-17, 9.1786093309558520e-20),
)).T
def planck(v: np.ndarray, T: float) -> np.ndarray:
return 2*h/c/c * v**3 / np.expm1(h*v/k/T)
guess = 2.5,
(T,), _ = curve_fit(
f=planck, xdata=freq, ydata=brightness_j, p0=guess, method='lm',
# bounds=(0.1, np.inf),
)
print('T =', T)
fig, ax = plt.subplots()
v = np.linspace(freq.min(), freq.max(), 500)
ax.scatter(freq, brightness_j, label='data')
ax.plot(v, planck(v, *guess), label='guess')
ax.plot(v, planck(v, T), label='fit')
ax.legend()
plt.show()
英文:
A lot of problems here, including that your variables were swapped, you're needlessly redefining physical constants, and your expression was highly numerically unstable. You need to use exp1m
instead:
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import h, c, k
from scipy.optimize import curve_fit
freq, brightness_erg, brightness_j = np.array((
(34473577711.372055, 7.0294715363905860e-16, 7.0294715363905860e-19),
(42896956937.695820, 1.0253178228238486e-15, 1.0253178228238486e-18),
(51322332225.447330, 1.3544045476166584e-15, 1.3544045476166584e-18),
(60344529880.182720, 1.6902073280174815e-15, 1.6902073280174815e-18),
(68767909106.506200, 2.0125779972022745e-15, 2.0125779972022743e-18),
(77780126454.101460, 2.3148004995630144e-15, 2.3148004995630145e-18),
(489996752265.52826, 3.2013198398211880e-16, 3.2013198398211880e-19),
(506039097962.67590, 2.5968748350997043e-16, 2.5968748350997040e-19),
(523273092332.36380, 2.0595903864583913e-16, 2.0595903864583912e-19),
(539918248580.78060, 1.7237876060575648e-16, 1.7237876060575649e-19),
(557158231134.75070, 1.3879848256567381e-16, 1.3879848256567383e-19),
(573803387383.16460, 1.0521820452559118e-16, 1.0521820452559118e-19),
(591049358121.42000, 9.1786093309558520e-17, 9.1786093309558520e-20),
)).T
def planck(v: np.ndarray, T: float) -> np.ndarray:
return 2*h/c/c * v**3 / np.expm1(h*v/k/T)
guess = 2.5,
(T,), _ = curve_fit(
f=planck, xdata=freq, ydata=brightness_j, p0=guess, method='lm',
# bounds=(0.1, np.inf),
)
print('T =', T)
fig, ax = plt.subplots()
v = np.linspace(freq.min(), freq.max(), 500)
ax.scatter(freq, brightness_j, label='data')
ax.plot(v, planck(v, *guess), label='guess')
ax.plot(v, planck(v, T), label='fit')
ax.legend()
plt.show()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论