英文:
What function would best fit the data I have from a galaxy?
问题
以下是代码的翻译部分：
我有以下一组数据：
surface_brightnesses_o2 = [12076.0616666451, 11850.730704516911, 10265.598145816548, 9120.859898168235, 7070.26133100111, 5636.138833975608, 3968.1608109082404, 2923.2839406153525, 1963.9315683870766, 1417.3534005331746, 953.9023540784231, 705.6331341427699, 494.19332394388607, 368.6833467905476, 266.41823769096874, 209.98748543636287, 162.17577134818487, 125.70474388251918, 99.72308185010249, 77.89696236284223, 53.44842864009773, 44.01192443651109, 35.52192383706094, 28.055033719366026]
surface_brightnesses_o3 = [24172.942124480545, 23257.99074788583, 19560.86193185194, 16867.86523112749, 12362.182457744273, 9447.974865736134, 6155.667579526176, 4233.309154367383, 2589.6992946467008, 1744.3756532539348, 1096.6861498588305, 768.600975237508, 512.7340397075068, 378.58271663510016, 268.4441550825379, 206.52758729119557, 155.45645416835472, 124.71693391104529, 97.34230151849876, 79.90134896492059, 63.519334039447266, 52.12382464229779, 41.91733978896593, 37.68365343589249, 31.54091147651983, 25.80764998552268, 22.808177293717083, 20.4718551088832, 16.05156984850126, 15.497358990115051, 15.42389243808505, 13.54177847744223]
radii_o2 = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0]
radii_o3 = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14, 14.5, 15, 15.5, 16]
surface_brightnesses_error_o2 = [109.89113552, 85.30012943, 80.8548183, 76.55283021, 66.49162753, 58.35388488, 49.4425817, 43.48019603, 36.48439283, 32.13758154, 28.57971998, 26.30618542, 24.27602806, 23.10048171, 22.01106869, 21.3172123, 20.77203895, 20.41962288, 20.12573286, 19.8928839, 19.84192745, 19.80754151, 19.6515864, 19.60323267]
surface_brightnesses_error_o3 = [155.47650023, 84.28555314, 74.17986129, 66.93258861, 54.67881726, 46.5099896, 36.86637245, 30.71396278, 25.45559327, 22.40018842, 19.83606727, 18.43327984, 16.94700871, 16.13059484, 15.55795461, 15.155422, 14.7707935, 14.59604581, 14.30144021, 14.13502224, 14.04555569, 13.9530354, 14.01473729, 14.13623735, 14.16959504, 14.1342218, 13.9836842, 13.87870645, 13.88701116, 13.91734777, 13.96048525, 13.98621865]
我试图绘制一个拟合曲线，使得y轴（surface brightnesses）是对数尺度，x轴（radii）是线性尺度。我还想在O2和O3的曲线中加入对应的误差。
我不想取surface brightness值的对数，只想按原样绘制数据并将y轴刻度设置为对数。然而，我找不到一个正确拟合数据的函数。
我会感激一些关于这里应该选择什么拟合以及如何编写代码的建议。
我尝试了拟合Sersic函数，这是一种用于研究星系表面亮度分布的亮度分布函数。
```python
fig, ax = plt.subplots(figsize=(10, 7))
# 定义Sersic函数
def sersic(r, I_e, R_e, n):
b_n = 1.9992 * n - 0.3271
return I_e * np.exp(-b_n * ((r/R_e)**(1/n) - 1))
# 对O2数据拟合模型
popt_o2, pcov_o2 = curve_fit(sersic, radii
<details>
<summary>英文:</summary>
I have the following set of data:
surface_brightnesses_o2 = [12076.0616666451, 11850.730704516911, 10265.598145816548, 9120.859898168235, 7070.26133100111, 5636.138833975608, 3968.1608109082404, 2923.2839406153525, 1963.9315683870766, 1417.3534005331746, 953.9023540784231, 705.6331341427699, 494.19332394388607, 368.6833467905476, 266.41823769096874, 209.98748543636287, 162.17577134818487, 125.70474388251918, 99.72308185010249, 77.89696236284223, 53.44842864009773, 44.01192443651109, 35.52192383706094, 28.055033719366026]
surface_brightnesses_o3 = [24172.942124480545, 23257.99074788583, 19560.86193185194, 16867.86523112749, 12362.182457744273, 9447.974865736134, 6155.667579526176, 4233.309154367383, 2589.6992946467008, 1744.3756532539348, 1096.6861498588305, 768.600975237508, 512.7340397075068, 378.58271663510016, 268.4441550825379, 206.52758729119557, 155.45645416835472, 124.71693391104529, 97.34230151849876, 79.90134896492059, 63.519334039447266, 52.12382464229779, 41.91733978896593, 37.68365343589249, 31.54091147651983, 25.80764998552268, 22.808177293717083, 20.4718551088832, 16.05156984850126, 15.497358990115051, 15.42389243808505, 13.54177847744223]
radii_o2 = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0]
radii_o3 = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14, 14.5, 15, 15.5, 16]
surface_brightnesses_error_o2 =
[109.89113552 85.30012943 80.8548183 76.55283021 66.49162753
58.35388488 49.4425817 43.48019603 36.48439283 32.13758154
28.57971998 26.30618542 24.27602806 23.10048171 22.01106869
21.3172123 20.77203895 20.41962288 20.12573286 19.8928839
19.84192745 19.80754151 19.6515864 19.60323267]
surface_brightnesses_error_o3 =
[155.47650023 84.28555314 74.17986129 66.93258861 54.67881726
46.5099896 36.86637245 30.71396278 25.45559327 22.40018842
19.83606727 18.43327984 16.94700871 16.13059484 15.55795461
15.155422 14.7707935 14.59604581 14.30144021 14.13502224
14.04555569 13.9530354 14.01473729 14.13623735 14.16959504
14.1342218 13.9836842 13.87870645 13.88701116 13.91734777
13.96048525 13.98621865]
I am trying to plot a fit such that the yscale (surface brightnesses) is log and the xscale (radii) is linear. I would also like to incorporate the errors for O2 and O3 in the corresponding plots for the surface brightnesses of O2 and O3.
I do not want to take log of the surface brightness values, I just want to plot the data as it is and set the yscale to log. However, I couldn't find a function that fits the data correctly.
I would appreciate some input on what would be a good fit here, and how to code it in.
I tried fitting a Sersic function, which is a brightness profile function used to study the surface brightness profiles of galaxies.
fig, ax = plt.subplots(figsize=(10, 7))
Define Sersic function
def sersic(r, I_e, R_e, n):
b_n = 1.9992*n - 0.3271
return I_e * np.exp(-b_n * ((r/R_e)**(1/n) - 1))
Fit the model to the O2 data
popt_o2, pcov_o2 = curve_fit(sersic, radii_o2, surface_brightnesses_o2, sigma=surface_brightnesses_error_o2, p0=[100000, 16, 2])
Fit the model to the O3 data
popt_o3, pcov_o3 = curve_fit(sersic, radii_o3, surface_brightnesses_o3, sigma=surface_brightnesses_error_o3, p0=[10000, 16, 2])
O2 data with error bars and fitted line
plt.errorbar(radii_o2, surface_brightnesses_o2, yerr=surface_brightnesses_error_o2, fmt='o', label='O2 data', capsize=4)
plt.plot(radii_o2, sersic(radii_o2, *popt_o2), 'r-', label='O2 fit: I_e=%5.3f, R_e=%5.3f, n=%5.3f' % tuple(popt_o2), color = 'blue')
O3 data with error bars and fitted line
plt.errorbar(radii_o3, surface_brightnesses_o3, yerr=surface_brightnesses_error_o3, fmt='o', label='O3 data', capsize=4)
plt.plot(radii_o3, sersic(radii_o3, *popt_o3), 'b-', label='O3 fit: I_e=%5.3f, R_e=%5.3f, n=%5.3f' % tuple(popt_o3), color = 'red')
plt.xlabel('Radii')
plt.ylabel('Surface Brightness')
plt.yscale('log')
plt.ylim(1, 30000) # Adjust the y-axis limits here
plt.title('Sersic Fit to Surface Brightness vs Radii for O2 and O3')
plt.legend()
plt.show()
And then I tried fitting a log-Gaussian plot:
Define the log-Gaussian function to fit to the data
def log_gaussian(x, amp, cen, wid):
return amp * np.exp(-(np.log(x) - cen)2 / wid2)
Initial guess for parameters (necessary for log-Gaussian)
popt_o2, pcov_o2 = curve_fit(power_law, radii_o2, surface_brightnesses_o2)
popt_o3, pcov_o3 = curve_fit(power_law, radii_o3, surface_brightnesses_o3
Fit the log-Gaussian model to the data
params_o2, _ = curve_fit(log_gaussian, radii_o2, surface_brightnesses_o2, p0_o2)
params_o3, _ = curve_fit(log_gaussian, radii_o3, surface_brightnesses_o3, p0_o3)
Generate points for the fitted log-Gaussian function
fit_o2 = power_law(radii_smooth_o2, *popt_o2)
fit_o3 = power_law(radii_smooth_o3, *popt_o3)
Create the plot
plt.figure(figsize=(10, 6))
Plot the original data
plt.errorbar(radii_o2, surface_brightnesses_o2, yerr=surface_brightnesses_error_o2, fmt='o', label='Data O2', capsize=4)
plt.errorbar(radii_o3, surface_brightnesses_o3, yerr=surface_brightnesses_error_o3, fmt='o', label='Data O3', capsize=4)
Plot the fitted log-Gaussian function
plt.plot(radii_fit, fit_o2, label='Fit O2', color = 'blue')
plt.plot(radii_fit, fit_o3, label='Fit O3', color = 'red')
Decorate the plot and set yscale to log
plt.xlabel('Radii')
plt.ylabel('Surface Brightnesses')
plt.title('Surface Brightnesses vs Radii')
plt.legend()
plt.yscale('log')
Show the plot
plt.show()
</details>
# 答案1
**得分**: 2
以下是代码的翻译部分：
1. 请使用不同的模型，并在这样做时执行对数拟合。
2. 在拟合期间，我认为您应该在y上应用对数，而不是x上应用。
有无限多个模型可供选择；哪些在科学上有效由您来决定。其中一个具有松散合理拟合的是具有线性衰减项的广义高斯模型；还有其他模型。
```python
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
surface_brightnesses_o2 = np.array([
12076.0616666451, 11850.730704516911, 10265.598145816548, 9120.859898168235, 7070.26133100111,
5636.138833975608, 3968.1608109082404, 2923.2839406153525, 1963.9315683870766, 1417.3534005331746,
953.9023540784231, 705.6331341427699, 494.19332394388607, 368.6833467905476, 266.41823769096874,
209.98748543636287, 162.17577134818487, 125.70474388251918, 99.72308185010249, 77.89696236284223,
53.44842864009773, 44.01192443651109, 35.52192383706094, 28.055033719366026
])
surface_brightnesses_o3 = np.array([
24172.942124480545, 23257.99074788583, 19560.86193185194, 16867.86523112749, 12362.182457744273,
9447.974865736134, 6155.667579526176, 4233.309154367383, 2589.6992946467008, 1744.3756532539348,
1096.6861498588305, 768.600975237508, 512.7340397075068, 378.58271663510016, 268.4441550825379,
206.52758729119557, 155.45645416835472, 124.71693391104529, 97.34230151849876, 79.90134896492059,
63.519334039447266, 52.12382464229779, 41.91733978896593, 37.68365343589249, 31.54091147651983,
25.80764998552268, 22.808177293717083, 20.4718551088832, 16.05156984850126, 15.497358990115051,
15.42389243808505, 13.54177847744223
])
radii_o2 = np.array([
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5,
10.0, 10.5, 11.0, 11.5, 12.0
])
radii_o3 = np.array([
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5,
10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14, 14.5, 15, 15.5, 16
])
surface_brightnesses_error_o2 = [
109.89113552, 85.30012943, 80.85481830, 76.55283021, 66.49162753,
58.353884880, 49.44258170, 43.48019603, 36.48439283, 32.13758154,
28.579719980, 26.30618542, 24.27602806, 23.10048171, 22.01106869,
21.317212300, 20.77203895, 20.41962288, 20.12573286, 19.89288390,
19.841927450, 19.80754151, 19.65158640, 19.60323267]
surface_brightnesses_error_o3 = [
155.47650023, 84.28555314, 74.17986129, 66.93258861, 54.67881726,
46.50998960, 36.86637245, 30.71396278, 25.45559327, 22.40018842,
19.83606727, 18.43327984, 16.94700871, 16.13059484, 15.55795461,
15.15542200, 14.77079350, 14.59604581, 14.30144021, 14.13502224,
14.04555569, 13.95303540, 14.01473729, 14.13623735, 14.16959504,
14.13422180, 13.98368420, 13.87870645, 13.88701116, 13.91734777,
13.96048525, 13.98621865]
def gaussian(x: np.ndarray, amp: float, cen: float, wid: float, pow: float, slope: float, off: float) -> np.ndarray:
return amp * np.exp(-np.abs((x - cen)/wid)**pow) + slope*x + off
def log_gaussian(x: np.ndarray, *params: float) -> np.ndarray:
return np.log(gaussian(x, *params))
ax: plt.Axes
fig, ax = plt
<details>
<summary>英文:</summary>
Use a different model, and when you do, perform a log-fit. You've applied your log on x when I believe you should apply it on y during fit.
There's an infinite number of models to choose from; which are scientifically valid is up to you to determine. One that has a loosely reasonable fit is a generalized Gaussian with a linear decay term; there are others.
```python
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
surface_brightnesses_o2 = np.array([
12076.0616666451, 11850.730704516911, 10265.598145816548, 9120.859898168235, 7070.26133100111,
5636.138833975608, 3968.1608109082404, 2923.2839406153525, 1963.9315683870766, 1417.3534005331746,
953.9023540784231, 705.6331341427699, 494.19332394388607, 368.6833467905476, 266.41823769096874,
209.98748543636287, 162.17577134818487, 125.70474388251918, 99.72308185010249, 77.89696236284223,
53.44842864009773, 44.01192443651109, 35.52192383706094, 28.055033719366026
])
surface_brightnesses_o3 = np.array([
24172.942124480545, 23257.99074788583, 19560.86193185194, 16867.86523112749, 12362.182457744273,
9447.974865736134, 6155.667579526176, 4233.309154367383, 2589.6992946467008, 1744.3756532539348,
1096.6861498588305, 768.600975237508, 512.7340397075068, 378.58271663510016, 268.4441550825379,
206.52758729119557, 155.45645416835472, 124.71693391104529, 97.34230151849876, 79.90134896492059,
63.519334039447266, 52.12382464229779, 41.91733978896593, 37.68365343589249, 31.54091147651983,
25.80764998552268, 22.808177293717083, 20.4718551088832, 16.05156984850126, 15.497358990115051,
15.42389243808505, 13.54177847744223
])
radii_o2 = np.array([
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5,
10.0, 10.5, 11.0, 11.5, 12.0
])
radii_o3 = np.array([
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5,
10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14, 14.5, 15, 15.5, 16
])
surface_brightnesses_error_o2 = [
109.89113552, 85.30012943, 80.85481830, 76.55283021, 66.49162753,
58.353884880, 49.44258170, 43.48019603, 36.48439283, 32.13758154,
28.579719980, 26.30618542, 24.27602806, 23.10048171, 22.01106869,
21.317212300, 20.77203895, 20.41962288, 20.12573286, 19.89288390,
19.841927450, 19.80754151, 19.65158640, 19.60323267]
surface_brightnesses_error_o3 = [
155.47650023, 84.28555314, 74.17986129, 66.93258861, 54.67881726,
46.50998960, 36.86637245, 30.71396278, 25.45559327, 22.40018842,
19.83606727, 18.43327984, 16.94700871, 16.13059484, 15.55795461,
15.15542200, 14.77079350, 14.59604581, 14.30144021, 14.13502224,
14.04555569, 13.95303540, 14.01473729, 14.13623735, 14.16959504,
14.13422180, 13.98368420, 13.87870645, 13.88701116, 13.91734777,
13.96048525, 13.98621865]
def gaussian(x: np.ndarray, amp: float, cen: float, wid: float, pow: float, slope: float, off: float) -> np.ndarray:
return amp * np.exp(-np.abs((x - cen)/wid)**pow) + slope*x + off
def log_gaussian(x: np.ndarray, *params: float) -> np.ndarray:
return np.log(gaussian(x, *params))
ax: plt.Axes
fig, ax = plt.subplots()
for title, brightness, radii, error, guess in (
(
'O2', surface_brightnesses_o2, radii_o2, surface_brightnesses_error_o2,
(1e4, 0, 1, 2, 0, 0),
),
(
'O3', surface_brightnesses_o3, radii_o3, surface_brightnesses_error_o3,
(1e4, 0, 1, 2, 0, 0),
),
):
ax.errorbar(radii, brightness, yerr=error, fmt='o', capsize=4, label=f'{title} data')
# ax.plot(radii, gaussian(radii, *guess), label=f'{title} guess')
fit, _ = curve_fit(
f=log_gaussian, xdata=radii, ydata=np.log(brightness), p0=guess,
bounds=(
( 1, -20, 0.01, 0.5, -1e6, -1e6),
(1e9, 20, 10.00, 10.0, 1e6, 1e6),
),
)
print(fit)
plot_radii = np.linspace(start=fit[1], stop=max(radii_o2.max(), radii_o3.max()), num=200)
ax.plot(plot_radii, gaussian(plot_radii, *fit), label=f'{title} fit')
plt.title('Gaussian Fit to Surface Brightness vs Radii for O2 and O3')
ax.set_xlabel('Radii')
ax.set_ylabel('Surface brightness')
ax.set_yscale('log')
ax.set_ylim(1, 1e5)
ax.legend()
plt.show()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论