在直方图上绘制正态分布

Question

我是 python 新手，在下面的代码中，我想绘制一条钟形曲线来显示数据如何遵循范数分布。我该怎么办？另外，任何人都可以回答为什么在显示 hist 时，我的值（x 轴）大于 100？我假设通过将 Randels 定义为 100，它不会显示任何高于它的东西。如果我没记错的话，x 轴代表我所在的“楼层”，y 轴代表有多少观察结果与该楼层匹配。顺便说一句，这是一个数据营项目。

"""
Let's say I roll a dice to determine if I go up or down a step in a building with
100 floors (1 step = 1 floor). If the dice is less than 2, I go down a step. If 
the dice is less than or equal to 5, I go up a step, and if the dice is equal to 6,
I go up x steps based on a random integer generator between 1 and 6. What is the probability
I will be higher than floor 60?
"""

import numpy as np
import matplotlib.pyplot as plt

# Set the seed
np.random.seed(123)

# Simulate random walk 
all_walks = []
for i in range(1000) :
    random_walk = [0]
    for x in range(100) :
        step = random_walk[-1]
        dice = np.random.randint(1,7)
        if dice <= 2:
            step = max(0, step - 1)
        elif dice <= 5:
            step = step + 1
        else:
            step = step + np.random.randint(1,7)
        if np.random.rand() <= 0.001 : # There's a 0.1% chance I fall and have to start at 0
            step = 0
        random_walk.append(step)
    all_walks.append(random_walk)

# Create and plot np_aw_t
np_aw_t = np.transpose(np.array(all_walks))

# Select last row from np_aw_t: ends
ends = np_aw_t[-1,:]

# Plot histogram of ends, display plot
plt.hist(ends,bins=10,edgecolor='k',alpha=0.65)
plt.style.use('fivethirtyeight')
plt.xlabel("Floor")
plt.ylabel("# of times in floor")
plt.show()

Answer 1

您可以使用scipy.stats.norm 获得正态分布。它的文档here。要将任何函数拟合到数据集，您可以使用scipy.optimize.curve_fit()，该文档为here。我的建议是这样的：

import scipy.stats as ss
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

#Making a figure with two y-axis (one for the hist, one for the pdf)
#An alternative would be to multiply the pdf by the sum of counts if you just want to show the fit.
fig, ax = plt.subplots(1,1)
twinx = ax.twinx()

rands = ss.norm.rvs(loc = 1, scale = 1, size = 1000)

#hist returns the bins and the value of each bin, plot to the y-axis ax
hist = ax.hist(rands)
vals, bins = hist[0], hist[1]

#calculating the center of each bin
bin_centers = [(bins[i] + bins[i+1])/2 for i in range(len(bins)-1)]

#finding the best fit coefficients, note vals/sum(vals) to get the probability in each bin instead of the count
coeff, cov = opt.curve_fit(ss.norm.pdf, bin_centers, vals/sum(vals), p0 = [0,1] )

#loc and scale are mean and standard deviation i believe
loc, scale = coeff

#x-values to plot the normal distribution curve
x = np.linspace(min(bins), max(bins), 100)

#Evaluating the pdf with the best fit mean and std
p = ss.norm.pdf(x, loc = loc, scale = scale)

#plot the pdf to the other axis and show
twinx.plot(x,p)
plt.show()

可能有更优雅的方法可以做到这一点，但如果您是 python 新手并且打算将它用于计算等，建议了解curve_fit 和scipy.stats。我不确定我是否理解“定义 Randels”的意思， hist 将绘制一个“标准”直方图，其中 x 轴上的 bin 和 y 轴上每个 bin 中的计数。当使用这些计数来拟合 pdf 时，我们可以将所有计数除以计数总数。

希望对你有帮助，如果有什么不清楚的地方就问吧:)

编辑：精简版

vals, bins,_ = ax.hist(my_histogram_data)
bin_centers = [(bins[i] + bins[i+1])/2 for i in range(len(bins)-1)]
coeff, cov = opt.curve_fit(ss.norm.pdf, bin_centers, vals/sum(vals), p0 = [0,1] )
x = np.linspace(min(bins), max(bins), 100)
p = ss.norm.pdf(x, loc = coeff[0], scale = coeff[1])
#p is now the fitted normal distribution