在python中将散点图变成直方图

Question

我需要根据另一个文件中的一些数据绘制直方图。

目前我有绘制散点图并拟合高斯的代码。

x 值是它正在读取的数据文件中相应行上的任何数字（在其他信息的前 12 行之后，即第 13 行是第一个事件），y 值是数字行数乘以一个值。

然后绘制并拟合散点图，但我需要能够将其绘制为直方图，并且能够更改 bin 宽度/数量（即将 bin 1、2、3 和 4 加在一起以获得 1/4整个垃圾箱的事件数量是事件数量的 4 倍——所以我猜想将数据中的多行加在一起），这就是我卡住的地方。

我将如何将其放入直方图并调整宽度/数字？

下面的代码，不知道如何使它漂亮。让我知道是否可以使它更易于阅读。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from numpy import exp, loadtxt, pi, sqrt, random, linspace
from lmfit import Model
import glob, os

## Define gaussian
def gaussian(x, amp, cen, wid):
    """1-d gaussian: gaussian(x, amp, cen, wid)"""
    return (amp / (sqrt(2*pi) * wid)) * exp(-(x-cen)**2 / (2*wid**2))

## Define constants
stderrThreshold = 10
minimumAmplitude = 0.1
approxcen = 780
MaestroT = 53

## Define paramaters
amps = []; ampserr = []; ts = []
folderToAnalyze = baseFolder + fileToRun + '\\'

## Generate the time array

for n in range(0, numfiles):
    
    ## Load text file
    x = np.linspace(0, 8191, 8192) 
    fullprefix = folderToAnalyze + prefix + str(n).zfill(3)
    y = loadtxt(fullprefix + ".Spe", skiprows= 12, max_rows = 8192) 

    ## Make figure
    fig, ax = plt.subplots(figsize=(15,8))
    fig.suptitle('Coincidence Detections', fontsize=20)
    plt.xlabel('Bins', fontsize=14)
    plt.ylabel('Counts', fontsize=14)

    ## Plot data
    ax.plot(x, y, 'bo')
    ax.set_xlim(600,1000)

    ## Fit data to Gaussian
    gmodel = Model(gaussian)
    result = gmodel.fit(y, x=x, amp=8, cen=approxcen, wid=1)

    ## Plot results and save figure
    ax.plot(x, result.best_fit, 'r-', label='best fit')
    ax.legend(loc='best')
    texttoplot = result.fit_report()
    ax.text(0.02, 0.5, texttoplot, transform=ax.transAxes)
    plt.close()
    fig.savefig(fullprefix + ".png", pad_inches='0.5')    

当前输出：散点图，确实显示了数据的预期分布和图（但是它们确实有一个糟糕的减少 chi^2，但一次有一个问题）

预期输出：相同数据的直方图，具有相同的分布和拟合，当每个事件被绘制为单独的 bin 时，希望可以将这些 bin 添加在一起以减少误差线

错误：不适用

数据：它基本上是超过 8192 行的标准分布。 1 个文件的完整数据为 here。还有原始的.Spe 文件、分散的plot 和完整版的code

2020-11-23 更新来自答案评论：

• 您好，我已经尝试实施了一段时间，但没有解决问题。我试图密切关注您的示例，但是我得到的直方图仍然具有 1 的 bin 宽度（即不加在一起）。我还在打印输出中获得了第二个空白图表，并且报告仅在 IDE 中打印输出（尽管我正在研究那个，并且估计我很快就会拥有它）。同样出于某种原因，它似乎在循环的 50 次迭代中的 3 次后停止。

• 这是当前状态下的code：

• 这是我得到的输出：

• 这是原始输出：

• 以防万一它有用，这是raw data。我似乎无法复制您的最后 2 个数字

• 理想的情况是能够将第 30 行的常数更改为所需的 bin 宽度，并在该情况下以该 bin 宽度运行。

Answer 1

• 在这种情况下，scatter plot 是一个直方图，除了点而不是条。
• .Spe 是每个事件的 bin 计数。
• x = np.linspace(0, 8191, 8192) 定义了 bins，bin 宽度为 1。
• 构建条形图而不是散点图
  • ax.bar(x, y) 而不是 ax.plot(x, y, 'bo')
• 由于现有数据，下图是分布非常广泛的直方图。
  • 取值范围从 321 到 1585
  • ax.set_xlim(300, 1800)

• 这些数据的好处是，可以很容易地根据x、bin 大小为 1 和y 重新创建原始分布，y 是每个x 的相应计数。
• np.repeat 可以创建一个包含重复元素的数组

import numpy
import matplotlib.pyplot as plt

# given x and y from the loop
# set the type as int
y = y.astype(int)
x = x.astype(int)

# create the data
data = np.repeat(x, y)

# determine the range of x
x_range = range(min(data), max(data)+1)

# determine the length of x
x_len = len(x_range)

# plot
fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(10, 10))

ax1.hist(data, bins=x_len)  # outliers are not plotted
ax2.boxplot(data, vert=False)
plt.show()

• 鉴于data，您现在可以执行任何需要的分析。
• SO: Fit gaussian to noisy data with lmfit
• LMFIT Docs
• Cross Validated 可能是深入研究模型的更好站点
• 所有误差计算参数均来自模型result。如果您从 np.histogram 为不同的 bin 宽度计算新的 x 和 y，这可能会影响错误。
  • approxcen = 780 也是result 的输入

# given x_len determine how many bins for a given bin width
width = 8
bins = int(np.round(x_len / width))

# determine new x and y for the histogram
y, x = np.histogram(data, bins=bins)

# Fit data to Gaussian
gmodel = Model(gaussian)
result = gmodel.fit(y, x=x[:-1], amp=8, cen=approxcen, wid=1)

# result
print(result.fit_report())

[out]:
[[Model]]
    Model(gaussian)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 314
    # data points      = 158
    # variables        = 3
    chi-square         = 397.702574
    reduced chi-square = 2.56582306
    Akaike info crit   = 151.851284
    Bayesian info crit = 161.039069
[[Variables]]
    amp:  1174.80608 +/- 37.1663147 (3.16%) (init = 8)
    cen:  775.535731 +/- 0.46232727 (0.06%) (init = 780)
    wid:  12.6563219 +/- 0.46232727 (3.65%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(amp, wid) =  0.577

# plot
plt.figure(figsize=(10, 6))
plt.bar(x[:-1], y)
plt.plot(x[:-1], result.best_fit, 'r-', label='best fit')

plt.figure(figsize=(20, 8))
plt.bar(x[:-1], y)
plt.xlim(700, 850)
plt.plot(x[:-1], result.best_fit, 'r-', label='best fit')
plt.grid()

• 从下一个代码块我们可以看出，错误与以下参数有关
  • stderrThreshold = 10
  • minimumAmplitude = 0.1
  • MaestroT = 53

    ## Append to list if error in amplitude and amplitude itself is within reasonable bounds
    if result.params['amp'].stderr < stderrThreshold and result.params['amp'] > minimumAmplitude:
        amps.append(result.params['amp'].value) 
        ampserr.append(result.params['amp'].stderr) 
        ts.append(MaestroT*n)