python - 如何使用直方图绘制状态密度，其曲线遵循python中每个峰的中点？

Question

我想绘制如上图（右下）所示的状态密度，其中直方图的一般形状显示在曲线中。该曲线遵循由峰的中点形成的一般形状。不幸的是，我只能像右上图那样绘制 DOS，其中曲线看起来像直方图本身。有没有办法用对应于直方图峰值中点的点来拟合曲线？

图片来自这些讲义。Eigenstates: [ 0.81 0.81 0.81 0.75 0.75 0.75 0.54 0.54 0.54 0.47 0.47 0.47 0.33 0.33 0.75 0.81 0.27 0.27 0.33 0.15 0.15 0.27 0.54 0.15 0.47 0.33 -0.06 -0.06 -0.06 0.27 0.15 -0.08 -0.06 -0.2 -0.2 -0.2 -1.03 - 0.08 -0.08 -0.27 -0.2 -1.03 -1.03 -1.03 -0.27 -0.27 -0.93 -0.93 -0.36 -0.36 -0.4 -0.4 -0.4 -0.27 -0.93 -0.36 -0.36 -0.4 -0.93 -0.93 -1.03 -1.03 - 0.93 -1.03 -1.03 -1.03 -1.03 -0.93 -0.08 -0.93 -1.03 -1.03 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 - 0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75]

我的代码：

bins = 100

plt.figure(figsize=[8,5])
plt.rcParams['axes.titlesize'] = 17
plt.rcParams.update({'font.size': 15})
dN = plt.hist(Eig_nn, bins)
plt.xlabel('Eigenvalues')
plt.ylabel('N')
plt.title('All eigenvalues')
plt.show()

# DOS

dE = (max(Eig_nn) - min(Eig_nn)) / bins
DOS = dN[0]/dE
E = np.linspace(-6,4,bins)

plt.figure(figsize=[8,5])
plt.rcParams['axes.titlesize'] = 17
plt.rcParams.update({'font.size': 15})
plt.plot(E, DOS)
plt.xlabel('E')
plt.ylabel('DOS')
plt.title('All eigenvalues')
plt.show()

我目前的尝试

标签： pythonmatplotlib

Answer 1

我不确定演示图是如何生成的。也许它只是用平滑曲线连接直方图的顶部？

它看起来有点像kde。这样的曲线可以通过 scipy's 生成gaussian_kde。它的第一个参数是数据点，第二个参数确定每个高斯的宽度。最佳宽度取决于数据以及您希望如何解释它们。对于此示例，大约 0.1 的值似乎足够了。目前还不清楚如何最好地将直方图与 kde 对齐，因此下面的代码使用辅助 y 轴绘制 kde。

一些实验代码：

from matplotlib import pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde

Eig_nn = np.array([...])

bins = 100
plt.rcParams.update({'font.size': 15})
fig, ax1 = plt.subplots(figsize=[8,5])
bin_values, bin_lims, _patches  = plt.hist(Eig_nn, bins, color='dodgerblue')

x = np.linspace(bin_lims[0]-0.2, bin_lims[-1]+0.2, 500) # create an x-axis
kde = gaussian_kde(Eig_nn, 0.1)
ax2 = ax1.twinx()  # create a secondary axis for the kde
ax2.plot(x, kde(x), color='crimson')

ax1.set_xlim(x[0], x[-1])  # set strict limits
ax1.set_xlabel('Eigenvalues')
ax1.set_ylabel('N')
ax2.set_ylabel('kde')
ax2.set_ylim(ymin=0) # put the zero of the secondary y-axis at the bottom
plt.title('All eigenvalues', fontsize=20)
plt.show()