使用scipy在xy平面定义一个monte carlo采样椭圆
Use scipy to define a monte carlo sampling ellipse in the xy plane
从 x--y 平面中的大量点开始，我想 select 这些点的一个子集，以便优先以已知主椭圆定义的方式和短轴。例如：
import numpy as np
npts = int(1e5)
lim = 3
x = np.random.uniform(-lim, lim, npts)
y = np.random.uniform(-lim, lim, npts)
major_axis = np.array((1, 1))
minor_axis = np.array((-0.25, 0.25))
以上两个向量定义了一个轴比为 4-1 的椭圆，长轴指向直线 y = x。因此，我试图写下一个 Monte Carlo 采样算法，如果 x--y 平面中的一个点位于输入长轴上（在本例中为 y=x 线），则它有更高的概率成为selected 在短轴上的一个点上（在本例中为 y = -x 线），其中概率增强因子仅由长轴与短轴之比确定（因此因子为 4 in这种情况）。
我一直在尝试使用scipy.stats.multivariate_normal的pdf方法来做到这一点，但我认为我一定是使用方法不当。我要解决的方法是通过将长轴和短轴视为特征方向来定义协方差矩阵，在每个点上使用 pdf 方法，对这些概率进行排序，然后 select 顶部 Nselect 这些概率。
from scipy.stats import multivariate_normal
cov = np.array((major_axis, minor_axis))
p = np.vstack((x, y)).T
prob_select = multivariate_normal.pdf(p, cov=cov)
idx_select = np.argsort(prob_select)
Nselect = len(x)/10
result_x = x[idx_select][-Nselect:]
result_y = y[idx_select][-Nselect:]
fig, ax = plt.subplots(1, 1)
__=ax.scatter(result_x, result_y, s=1)
xlim = ax.set_xlim(-3, 3)
ylim = ax.set_ylim(-3, 3)
上图表明我的算法有问题，因为这个椭圆的长轴不在 y=x 线上。我怀疑协方差矩阵没有正确定义，但是当我将相同的协方差矩阵与 rvs 方法一起使用时，我得到了预期的分布：
correct_result = multivariate_normal.rvs(size=Nselect, cov=cov)
fig, ax = plt.subplots(1, 1)
__=ax.scatter(correct_result[:, 0], correct_result[:, 1], s=1)
xlim = ax.set_xlim(-3, 3)
ylim = ax.set_ylim(-3, 3)
我在使用 multivariate_normal.pdf 或协方差矩阵定义时是否犯了一个简单的错误？如果算法在其他方面存在缺陷，是否有更简单的方法来定义从椭圆的 major/minor 轴开始的 selection 函数？
为什么结果不一致
此处的协方差矩阵格式错误，您无法从结果行为中做出任何推论。 rvs 方法在这种情况下给出不同结果的事实仅反映了 rvs 和 pdf 函数以不同方式预处理它们的参数这一事实。而 rvs 基本上将其参数直接传递给 numpy.multivariate_normal...
# https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/_multivariate.py#L405
dim, mean, cov = _process_parameters(None, mean, cov)
out = np.random.multivariate_normal(mean, cov, size)
return _squeeze_output(out)
pdf 将协方差矩阵传递给计算伪逆的函数：
# https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/_multivariate.py#L378
dim, mean, cov = _process_parameters(None, mean, cov)
x = _process_quantiles(x, dim)
prec_U, log_det_cov = _psd_pinv_decomposed_log_pdet(cov)
out = np.exp(self._logpdf(x, mean, prec_U, log_det_cov))
return _squeeze_output(out)
如果协方差矩阵格式正确，这些只能保证给出一致的结果。
建立一个结构良好的矩阵
协方差矩阵只是每个对应维度对的协方差，因此根据定义它是对称的。
文档重申了这一点：
cov : 2-D array_like, of shape (N, N)
Covariance matrix of the distribution. It must be symmetric and positive-semidefinite for proper sampling.
假设你想得到一个长轴和短轴的协方差矩阵，你真正想要的是求解一个反向特征向量问题！耶！我希望我们有 mathjax...
我们需要一个对称矩阵C = [[a, b], [b, a]]，使得[1, 1]和[1, -1]是特征向量，我们还希望特征值的比例是四比一。这意味着 C * [1, 1] = [4, 4] 和 C * [1, -1] = [1, -1]。选择 1 作为我们的次要索引特征值和 4 作为我们的主要索引特征值，并使用变量手动进行矩阵乘法，我们得到 a + b = 4 和 a - b = 1。所以A是2.5，b是1.5，C是[[2.5, 1.5], [1.5, 2.5]].
我们还可以使用矩阵方程来找到更直接的解。如果 E 是特征向量 [[1, 1], [1, -1]] 的矩阵并且 lambda 是特征值 [[4, 0], [0, 1]] 的对角矩阵，那么我们正在寻找矩阵 X 使得：
X @ E = E @ lambda
其中 @ 表示矩阵乘法（如 Python 3.5+）。
也就是说
X = E @ lambda @ E ^ -1
在numpy中是
>>> E = numpy.array([[1, 1], [1, -1]])
>>> lambda_ = numpy.array([[4, 0], [0, 1]])
>>> E @ lambda_ @ numpy.linalg.pinv(E)
array([[ 2.5, 1.5],
[ 1.5, 2.5]])
在您的代码中将其用作 cov 可得到以下结果：
从 x--y 平面中的大量点开始，我想 select 这些点的一个子集，以便优先以已知主椭圆定义的方式和短轴。例如：
import numpy as np
npts = int(1e5)
lim = 3
x = np.random.uniform(-lim, lim, npts)
y = np.random.uniform(-lim, lim, npts)
major_axis = np.array((1, 1))
minor_axis = np.array((-0.25, 0.25))
以上两个向量定义了一个轴比为 4-1 的椭圆，长轴指向直线 y = x。因此，我试图写下一个 Monte Carlo 采样算法，如果 x--y 平面中的一个点位于输入长轴上（在本例中为 y=x 线），则它有更高的概率成为selected 在短轴上的一个点上（在本例中为 y = -x 线），其中概率增强因子仅由长轴与短轴之比确定（因此因子为 4 in这种情况）。
我一直在尝试使用scipy.stats.multivariate_normal的pdf方法来做到这一点，但我认为我一定是使用方法不当。我要解决的方法是通过将长轴和短轴视为特征方向来定义协方差矩阵，在每个点上使用 pdf 方法，对这些概率进行排序，然后 select 顶部 Nselect 这些概率。
from scipy.stats import multivariate_normal
cov = np.array((major_axis, minor_axis))
p = np.vstack((x, y)).T
prob_select = multivariate_normal.pdf(p, cov=cov)
idx_select = np.argsort(prob_select)
Nselect = len(x)/10
result_x = x[idx_select][-Nselect:]
result_y = y[idx_select][-Nselect:]
fig, ax = plt.subplots(1, 1)
__=ax.scatter(result_x, result_y, s=1)
xlim = ax.set_xlim(-3, 3)
ylim = ax.set_ylim(-3, 3)
上图表明我的算法有问题，因为这个椭圆的长轴不在 y=x 线上。我怀疑协方差矩阵没有正确定义，但是当我将相同的协方差矩阵与 rvs 方法一起使用时，我得到了预期的分布：
correct_result = multivariate_normal.rvs(size=Nselect, cov=cov)
fig, ax = plt.subplots(1, 1)
__=ax.scatter(correct_result[:, 0], correct_result[:, 1], s=1)
xlim = ax.set_xlim(-3, 3)
ylim = ax.set_ylim(-3, 3)
我在使用 multivariate_normal.pdf 或协方差矩阵定义时是否犯了一个简单的错误？如果算法在其他方面存在缺陷，是否有更简单的方法来定义从椭圆的 major/minor 轴开始的 selection 函数？
为什么结果不一致
此处的协方差矩阵格式错误，您无法从结果行为中做出任何推论。 rvs 方法在这种情况下给出不同结果的事实仅反映了 rvs 和 pdf 函数以不同方式预处理它们的参数这一事实。而 rvs 基本上将其参数直接传递给 numpy.multivariate_normal...
# https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/_multivariate.py#L405
dim, mean, cov = _process_parameters(None, mean, cov)
out = np.random.multivariate_normal(mean, cov, size)
return _squeeze_output(out)
pdf 将协方差矩阵传递给计算伪逆的函数：
# https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/_multivariate.py#L378
dim, mean, cov = _process_parameters(None, mean, cov)
x = _process_quantiles(x, dim)
prec_U, log_det_cov = _psd_pinv_decomposed_log_pdet(cov)
out = np.exp(self._logpdf(x, mean, prec_U, log_det_cov))
return _squeeze_output(out)
如果协方差矩阵格式正确，这些只能保证给出一致的结果。
建立一个结构良好的矩阵
协方差矩阵只是每个对应维度对的协方差，因此根据定义它是对称的。
文档重申了这一点：
cov : 2-D array_like, of shape (N, N)
Covariance matrix of the distribution. It must be symmetric and positive-semidefinite for proper sampling.
假设你想得到一个长轴和短轴的协方差矩阵，你真正想要的是求解一个反向特征向量问题！耶！我希望我们有 mathjax...
我们需要一个对称矩阵C = [[a, b], [b, a]]，使得[1, 1]和[1, -1]是特征向量，我们还希望特征值的比例是四比一。这意味着 C * [1, 1] = [4, 4] 和 C * [1, -1] = [1, -1]。选择 1 作为我们的次要索引特征值和 4 作为我们的主要索引特征值，并使用变量手动进行矩阵乘法，我们得到 a + b = 4 和 a - b = 1。所以A是2.5，b是1.5，C是[[2.5, 1.5], [1.5, 2.5]].
我们还可以使用矩阵方程来找到更直接的解。如果 E 是特征向量 [[1, 1], [1, -1]] 的矩阵并且 lambda 是特征值 [[4, 0], [0, 1]] 的对角矩阵，那么我们正在寻找矩阵 X 使得：
X @ E = E @ lambda
其中 @ 表示矩阵乘法（如 Python 3.5+）。
也就是说
X = E @ lambda @ E ^ -1
在numpy中是
>>> E = numpy.array([[1, 1], [1, -1]])
>>> lambda_ = numpy.array([[4, 0], [0, 1]])
>>> E @ lambda_ @ numpy.linalg.pinv(E)
array([[ 2.5, 1.5],
[ 1.5, 2.5]])
在您的代码中将其用作 cov 可得到以下结果：