使用辛普森规则对正态分布进行积分答案

【问题标题】：Integration of normal distribution using Simpson's rule使用辛普森规则对正态分布进行积分
【发布时间】：2013-12-21 04:56:05
【问题描述】：

我正在尝试使用 scipy.integrate.simps 函数执行简单的集成，但我无法弄清楚它显示的结果。

这是一个 MWE：

import numpy as np
from scipy.integrate import simps

# Same normal function used by np.random.normal
def norm_func(x, mu, sigma):
    y = 1/(sigma*np.sqrt(2*np.pi))*np.exp(-(x-mu)**2/(2*sigma**2))
    return y

# Generate some random points from the normal distribution.
a = np.random.normal(1., 0.1, 1000)

# Integrate the evaluated values of these points.
print simps(norm_func(a, 1., 0.1), a)

我希望由于我从normal distribution 中抽取随机数，然后将它们的评估整合到等效的正态分布中，我应该得到整合所述正态分布which is 1 的结果（或非常接近它） .

我发现结果似乎随着a 的样本大小而变化。更糟糕的是，如果我在 a = np.random.normal(1., 0.1, 10000) 中设置 10000 的值，则集成返回一个负值。

我在这里做错了什么？

【问题讨论】：

标签： python numpy integration

【解决方案1】：

使用您的样本，首先对a 进行排序，因为它应该是一个要采样的点数组，并且它希望它们是为了构建近似值。辛普森法则使用

因此它将从您的数组中获取x 的值并评估该函数。如果它们的顺序是随机的，你会发现上面的公式没有什么意义，因为它会从域上的一个随机点到另一个点进行积分。将其视为x 可能会更好，所以我将使用该变量名：

x = np.random.normal(1., 0.1, 1000)
x.sort() # sorts in place
print simps(norm_func(x, 1., 0.1), x)
#0.999914876748

这也适用于我：

s = np.sort(np.random.normal(1., 0.1, 10000))
print simps(norm_func(s, 1., 0.1), s)
#0.999943377731

【讨论】：

我确实尝试对列表进行排序，但我使用了sorted(a) 而不是a.sort()，但它不起作用。知道为什么会这样吗？
也许你没有保存sorted的输出？ a.sort() 对 a 就地排序，所以 a.sort() 等价于 a = np.sort(a) 或 a = np.array(sorted(a))。请参阅我的编辑以了解如何做到这一点。
不，我确实保存了输出。你能不能自己试一试，看看是不是只有我一个人？这是我对列表进行排序的操作：a1 = sorted(a) / print simps(norm_func(a, 1., 0.1), a1).
从头开始，我刚刚意识到我只将最后一个 a 更改为 a1 而不是用于调用普通函数的那个。

【解决方案2】：

以下是使用simps 集成标准普通pdf 的方法：

In [37]: a = np.linspace(-20, 20)

In [38]: print simps(norm_func(a, 0., 1.), a)
1.0

【讨论】：