Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

huangapple
go评论93阅读模式
英文:

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

问题

协方差未在SciPy的Curvefit中估算

这是我的数据集：

  1. frequency (Hz)    brightness (ergs/s/cm^2/sr/Hz)    brightness (J/s/m^2/sr/Hz)
  2. float64    float64    float64
  3. 34473577711.372055    7.029471536390586e-16    7.029471536390586e-19
  4. 42896956937.69582    1.0253178228238486e-15    1.0253178228238486e-18
  5. 51322332225.44733    1.3544045476166584e-15    1.3544045476166584e-18
  6. 60344529880.18272    1.6902073280174815e-15    1.6902073280174815e-18
  7. 68767909106.5062    2.0125779972022745e-15    2.0125779972022743e-18
  8. 77780126454.10146    2.3148004995630144e-15    2.3148004995630145e-18
  9. ...    ...    ...
  10. 489996752265.52826    3.201319839821188e-16    3.201319839821188e-19
  11. 506039097962.6759    2.5968748350997043e-16    2.596874835099704e-19
  12. 523273092332.3638    2.0595903864583913e-16    2.0595903864583912e-19
  13. 539918248580.7806    1.7237876060575648e-16    1.7237876060575649e-19
  14. 557158231134.7507    1.3879848256567381e-16    1.3879848256567383e-19
  15. 573803387383.1646    1.0521820452559118e-16    1.0521820452559118e-19
  16. 591049358121.42    9.178609330955852e-17    9.178609330955852e-20

我尝试使用CurveFit将其拟合到普朗克的辐射定律：

  1. import numpy as np
  2. from scipy.optimize import curve_fit
  4. h = 6.626 * 10e-34
  5. c = 3 * 10e8
  6. k = 1.38 * 10e-23
  7. const1 = 2 * h / (c**2)
  8. const2 = h / k
  10. def planck(x, v):
  11.     return const1 * (v**3) * (1 / ((np.exp(const2 * v / x)) - 1))
  13. popt, pcov = curve_fit(planck, cmb['frequency (Hz)'], cmb['brightness (J/s/m^2/sr/Hz)'])
  14. print(popt, pcov)

警告：

  1. /tmp/ipykernel_2500/4072287013.py:11: RuntimeWarning: divide by zero encountered in divide
  2.   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:

  1. frequency (Hz)	brightness (ergs/s/cm^2/sr/Hz)	brightness (J/s/m^2/sr/Hz)
  2. float64	float64	float64
  3. 34473577711.372055	7.029471536390586e-16	7.029471536390586e-19
  4. 42896956937.69582	1.0253178228238486e-15	1.0253178228238486e-18
  5. 51322332225.44733	1.3544045476166584e-15	1.3544045476166584e-18
  6. 60344529880.18272	1.6902073280174815e-15	1.6902073280174815e-18
  7. 68767909106.5062	2.0125779972022745e-15	2.0125779972022743e-18
  8. 77780126454.10146	2.3148004995630144e-15	2.3148004995630145e-18
  9. ...	...	...
  10. 489996752265.52826	3.201319839821188e-16	3.201319839821188e-19
  11. 506039097962.6759	2.5968748350997043e-16	2.596874835099704e-19
  12. 523273092332.3638	2.0595903864583913e-16	2.0595903864583912e-19
  13. 539918248580.7806	1.7237876060575648e-16	1.7237876060575649e-19
  14. 557158231134.7507	1.3879848256567381e-16	1.3879848256567383e-19
  15. 573803387383.1646	1.0521820452559118e-16	1.0521820452559118e-19
  16. 591049358121.42	9.178609330955852e-17	9.178609330955852e-20

I tried to use CurveFit to fit this to Planck's Radiation Law:

  1. import numpy as np
  2. from scipy.optimize import curve_fit
  4. h=6.626*10e-34
  5. c=3*10e8
  6. k=1.38*10e-23
  7. const1=2*h/(c**2)
  8. const2=h/k
  10. def planck(x,v):
  11.     return const1*(v**3)*(1/((np.exp(const2*v/x))-1))
  13. popt,pcov= curve_fit(planck, cmb['frequency (Hz)'],cmb['brightness (J/s/m^2/sr/Hz)'])
  14. print(popt, pcov)

Warning:

  1. /tmp/ipykernel_2500/4072287013.py:11: RuntimeWarning: divide by zero encountered in divide
  2.   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代替：

  1. import matplotlib.pyplot as plt
  2. import numpy as np
  3. from scipy.constants import h, c, k
  4. from scipy.optimize import curve_fit
  6. freq, brightness_erg, brightness_j = np.array(((
  7.     (34473577711.372055, 7.0294715363905860e-16, 7.0294715363905860e-19),
  8.     (42896956937.695820, 1.0253178228238486e-15, 1.0253178228238486e-18),
  9.     (51322332225.447330, 1.3544045476166584e-15, 1.3544045476166584e-18),
  10.     (60344529880.182720, 1.6902073280174815e-15, 1.6902073280174815e-18),
  11.     (68767909106.506200, 2.0125779972022745e-15, 2.0125779972022743e-18),
  12.     (77780126454.101460, 2.3148004995630144e-15, 2.3148004995630145e-18),
  13.     (489996752265.52826, 3.2013198398211880e-16, 3.2013198398211880e-19),
  14.     (506039097962.67590, 2.5968748350997043e-16, 2.5968748350997040e-19),
  15.     (523273092332.36380, 2.0595903864583913e-16, 2.0595903864583912e-19),
  16.     (539918248580.78060, 1.7237876060575648e-16, 1.7237876060575649e-19),
  17.     (557158231134.75070, 1.3879848256567381e-16, 1.3879848256567383e-19),
  18.     (573803387383.16460, 1.0521820452559118e-16, 1.0521820452559118e-19),
  19.     (591049358121.42000, 9.1786093309558520e-17, 9.1786093309558520e-20),
  20. )).T
  22. def planck(v: np.ndarray, T: float) -> np.ndarray:
  23.     return 2*h/c/c * v**3 / np.expm1(h*v/k/T)
  25. guess = 2.5,
  27. (T,), _ = curve_fit(
  28.     f=planck, xdata=freq, ydata=brightness_j, p0=guess, method='lm',
  29.     # bounds=(0.1, np.inf),
  30. )
  31. print('T =', T)
  33. fig, ax = plt.subplots()
  34. v = np.linspace(freq.min(), freq.max(), 500)
  35. ax.scatter(freq, brightness_j, label='data')
  36. ax.plot(v, planck(v, *guess), label='guess')
  37. ax.plot(v, planck(v, T), label='fit')
  38. ax.legend()
  39. plt.show()

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

英文:

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:

  1. import matplotlib.pyplot as plt
  2. import numpy as np
  3. from scipy.constants import h, c, k
  4. from scipy.optimize import curve_fit
  6. freq, brightness_erg, brightness_j = np.array((
  7.     (34473577711.372055, 7.0294715363905860e-16, 7.0294715363905860e-19),
  8.     (42896956937.695820, 1.0253178228238486e-15, 1.0253178228238486e-18),
  9.     (51322332225.447330, 1.3544045476166584e-15, 1.3544045476166584e-18),
  10.     (60344529880.182720, 1.6902073280174815e-15, 1.6902073280174815e-18),
  11.     (68767909106.506200, 2.0125779972022745e-15, 2.0125779972022743e-18),
  12.     (77780126454.101460, 2.3148004995630144e-15, 2.3148004995630145e-18),
  13.     (489996752265.52826, 3.2013198398211880e-16, 3.2013198398211880e-19),
  14.     (506039097962.67590, 2.5968748350997043e-16, 2.5968748350997040e-19),
  15.     (523273092332.36380, 2.0595903864583913e-16, 2.0595903864583912e-19),
  16.     (539918248580.78060, 1.7237876060575648e-16, 1.7237876060575649e-19),
  17.     (557158231134.75070, 1.3879848256567381e-16, 1.3879848256567383e-19),
  18.     (573803387383.16460, 1.0521820452559118e-16, 1.0521820452559118e-19),
  19.     (591049358121.42000, 9.1786093309558520e-17, 9.1786093309558520e-20),
  20. )).T
  23. def planck(v: np.ndarray, T: float) -> np.ndarray:
  24.     return 2*h/c/c * v**3 / np.expm1(h*v/k/T)
  26. guess = 2.5,
  28. (T,), _ = curve_fit(
  29.     f=planck, xdata=freq, ydata=brightness_j, p0=guess, method='lm',
  30.     # bounds=(0.1, np.inf),
  31. )
  32. print('T =', T)
  34. fig, ax = plt.subplots()
  35. v = np.linspace(freq.min(), freq.max(), 500)
  36. ax.scatter(freq, brightness_j, label='data')
  37. ax.plot(v, planck(v, *guess), label='guess')
  38. ax.plot(v, planck(v, T), label='fit')
  39. ax.legend()
  40. plt.show()

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

huangapple
  • 本文由 发表于 2023年6月1日 21:28:27
  • 转载请务必保留本文链接：https://go.coder-hub.com/76382415.html
  • astronomy
  • python
  • scientific-computing
  • scipy
  • scipy-optimize
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定