在 python 中将散点图转换为直方图
Turning a scatter plot into a histogram in python
我需要根据另一个文件中的一些数据绘制直方图。
目前我有绘制散点图并适合高斯分布的代码。
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 加在一起以减少误差线
错误：N/A
数据：基本上是8192行的标准分布。 1 个文件的完整数据是 here. Also the original .Spe file, the scatter plot and the full version of the code
2020-11-23 回复评论更新：
您好，我已经尝试实现这一点一段时间了，但没有成功。我试图密切关注您的示例，但是我得到的直方图仍然具有 1 的 bin 宽度（即不加在一起）。我还在打印输出中得到第二张空白图表，并且报告仅在 IDE 中打印输出（尽管我正在处理那个，并且我认为我很快就会得到它）。同样出于某种原因，它似乎在 50 次循环迭代中的 3 次后停止。
这是 code 当前状态：
这是我得到的输出：
这是原始输出：
为了以防万一，这是 raw data。我似乎无法复制你的最后 2 个数字
理想的情况是能够将第 30 行上的常量更改为任何所需的 bin 宽度，并在该场合使用该 bin 宽度 运行。
- 在这种情况下，
scatter plot 是直方图，只是用点代替了条。
.Spe 是每个事件的 bin 计数。x = np.linspace(0, 8191, 8192)定义bin，bin宽度为1。- 构建条形图而不是散点图
ax.bar(x, y) 而不是 ax.plot(x, y, 'bo')
- 由于已有数据，下图是一个分布很广的直方图。
- 有从 321 到 1585 的值
ax.set_xlim(300, 1800)
- 此数据的好处是，可以很容易地重新创建基于
x 的原始分布，bin 大小为 1，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，您现在可以执行任何需要的分析。
- 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 = 10minimumAmplitude = 0.1MaestroT = 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)
我需要根据另一个文件中的一些数据绘制直方图。
目前我有绘制散点图并适合高斯分布的代码。
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 加在一起以减少误差线
错误：N/A
数据：基本上是8192行的标准分布。 1 个文件的完整数据是 here. Also the original .Spe file, the scatter plot and the full version of the code
2020-11-23 回复评论更新：
您好，我已经尝试实现这一点一段时间了，但没有成功。我试图密切关注您的示例，但是我得到的直方图仍然具有 1 的 bin 宽度（即不加在一起）。我还在打印输出中得到第二张空白图表，并且报告仅在 IDE 中打印输出（尽管我正在处理那个，并且我认为我很快就会得到它）。同样出于某种原因，它似乎在 50 次循环迭代中的 3 次后停止。
这是 code 当前状态：
这是我得到的输出：
这是原始输出：
为了以防万一，这是 raw data。我似乎无法复制你的最后 2 个数字
理想的情况是能够将第 30 行上的常量更改为任何所需的 bin 宽度，并在该场合使用该 bin 宽度 运行。
- 在这种情况下，
scatter plot是直方图，只是用点代替了条。 .Spe是每个事件的 bin 计数。x = np.linspace(0, 8191, 8192)定义bin，bin宽度为1。- 构建条形图而不是散点图
ax.bar(x, y)而不是ax.plot(x, y, 'bo')
- 由于已有数据，下图是一个分布很广的直方图。
- 有从 321 到 1585 的值
ax.set_xlim(300, 1800)
- 此数据的好处是，可以很容易地重新创建基于
x的原始分布，bin 大小为 1，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，您现在可以执行任何需要的分析。 - 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 = 10minimumAmplitude = 0.1MaestroT = 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)