我有以下一组数据:
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]
我试图绘制一个拟合图,使yscale(表面亮度)为log,xscale(半径)为线性。我还想将O2和O3的误差合并到相应的O2和O3的表面亮度图中。
我不想记录表面亮度值,我只想按原样绘制数据,并将yscale设置为log。但是,我找不到一个函数,可以正确地拟合数据。
我将感谢一些输入什么将是一个很好的适合在这里,以及如何编码它。
我试着拟合一个Sersic函数,这是一个用于研究星系表面亮度分布的亮度分布函数。
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()
然后我尝试拟合对数高斯图:
# 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 / wid**2)
# 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()
1条答案
按热度按时间yyhrrdl81#
使用不同的模型,然后执行对数拟合。你已经在x上应用了你的log,而我认为你应该在fit期间在y上应用它。
有无数的模型可供选择;哪些是科学上有效的由你来决定一个有一个松散合理的拟合是一个广义高斯与线性衰减项;还有其他人