加快正态分布概率质量分配

Question

我们有 N 个用户，平均为 P。每个用户的点数，其中每个点是 0 到 1 之间的单个值。我们需要使用已知密度为 0.05 的正态分布来分配每个点的质量，因为这些点具有一些不确定性。此外，我们需要将质量包裹在 0 和 1 周围，例如0.95 处的点也将分配 0 左右的质量。我在下面提供了一个工作示例，它将正态分布分箱到 D=50 箱中。该示例使用 Python 类型模块，但如果您愿意，可以忽略它。

from typing import List, Any
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt

D = 50
BINS: List[float] = np.linspace(0, 1, D + 1).tolist()


def probability_mass(distribution: Any, x0: float, x1: float) -> float:
    """
    Computes the area under the distribution, wrapping at 1.
    The wrapping is done by adding the PDF at +- 1.
    """
    assert x1 > x0
    return (
        (distribution.cdf(x1) - distribution.cdf(x0))
        + (distribution.cdf(x1 + 1) - distribution.cdf(x0 + 1))
        + (distribution.cdf(x1 - 1) - distribution.cdf(x0 - 1))
    )


def point_density(x: float) -> List[float]:
    distribution: Any = scipy.stats.norm(loc=x, scale=0.05)
    density: List[float] = []
    for i in range(D):
        density.append(probability_mass(distribution, BINS[i], BINS[i + 1]))
    return density


def user_density(points: List[float]) -> Any:

    # Find the density of each point
    density: Any = np.array([point_density(p) for p in points])

    # Combine points and normalize
    combined = density.sum(axis=0)
    return combined / combined.sum()


if __name__ == "__main__":

    # Example for one user
    data: List[float] = [.05, .3, .5, .5]
    density = user_density(data)

    # Example for multiple users (N = 2)
    print([user_density(x) for x in [[.3, .5], [.7, .7, .7, .9]]])

    ### NB: THE REMAINING CODE IS FOR ILLUSTRATION ONLY!
    ### NB: THE IMPORTANT THING IS TO COMPUTE THE DENSITY FAST!
    middle: List[float] = []
    for i in range(D):
        middle.append((BINS[i] + BINS[i + 1]) / 2)
    plt.bar(x=middle, height=density, width=1.0 / D + 0.001)
    plt.xlim(0, 1)
    plt.xlabel("x")
    plt.ylabel("Density")
    plt.show()

在此示例中为 N=1、D=50、P=4。但是，我们希望将这种方法扩展到 N=10000 和 P=100，同时尽可能快。我不清楚我们如何向量化这种方法。我们如何才能最好地加快速度？

编辑

更快的解决方案可能会产生略微不同的结果。例如，它可以近似正态分布，而不是使用精确的正态分布。

EDIT2

我们只关心使用user_density() 函数计算density。该情节仅用于帮助解释该方法。我们不关心情节本身:)

EDIT3

请注意，P 是平均值。每个用户的积分。一些用户可能拥有更多，而一些用户可能拥有更少。如果有帮助，您可以假设我们可以扔掉积分，这样所有用户的最大积分为 2 * P。只要解决方案可以为每个用户处理灵活的点数，就可以在进行基准测试时忽略这部分。

Answer 1

这将是我的矢量化方法：

data = np.array([0.05, 0.3, 0.5, 0.5])

np.random.seed(31415)
# random noise
randoms = np.random.normal(0,1,(len(data), int(1e5))) * 0.05

# samples with noise
samples = data[:,None] + randoms

# wrap [0,1]
samples = (samples % 1).ravel()


# histogram
hist, bins, patches = plt.hist(samples, bins=BINS, density=True)

输出：

Answer 2

我能够将时间从 100 个数据点的每个样本大约 4 秒减少到每个样本大约 1 毫秒。

在我看来，您花费了大量时间来模拟大量正态分布。由于无论如何您都在处理非常大的样本量，因此您不妨只使用标准正态分布值，因为无论如何它都会取平均值。

我重新创建了您的方法（BaseMethod 类），然后创建了一个优化类（OptimizedMethod 类），并使用 timeit 装饰器对它们进行了评估。我的方法的主要区别在于以下行：

    # Generate a standardized set of values to add to each sample to simulate normal distribution
    self.norm_vals = np.array([norm.ppf(x / norm_val_n) * 0.05 for x in range(1, norm_val_n, 1)])

这会根据逆正态累积分布函数创建一组通用数据点，我们可以将其添加到每个数据点以模拟该点周围的正态分布。然后我们只是将数据重塑为用户样本并在样本上运行 np.histogram。

import numpy as np
import scipy.stats
from scipy.stats import norm
import time

# timeit decorator for evaluating performance
def timeit(method):
    def timed(*args, **kw):
        ts = time.time()
        result = method(*args, **kw)
        te = time.time()
        print('%r  %2.2f ms' % (method.__name__, (te - ts) * 1000 ))
        return result
    return timed

# Define Variables
N = 10000
D = 50
P = 100

# Generate sample data
np.random.seed(0)
data = np.random.rand(N, P)

# Run OP's method for comparison
class BaseMethod:
    def __init__(self, d=50):
        self.d = d
        self.bins = np.linspace(0, 1, d + 1).tolist()

    def probability_mass(self, distribution, x0, x1):
        """
        Computes the area under the distribution, wrapping at 1.
        The wrapping is done by adding the PDF at +- 1.
        """
        assert x1 > x0
        return (
            (distribution.cdf(x1) - distribution.cdf(x0))
            + (distribution.cdf(x1 + 1) - distribution.cdf(x0 + 1))
            + (distribution.cdf(x1 - 1) - distribution.cdf(x0 - 1))
        )

    def point_density(self, x):
        distribution = scipy.stats.norm(loc=x, scale=0.05)
        density = []
        for i in range(self.d):
            density.append(self.probability_mass(distribution, self.bins[i], self.bins[i + 1]))
        return density

    @timeit
    def base_user_density(self, data):
        n = data.shape[0]
        density = np.empty((n, self.d))
        for i in range(data.shape[0]):
            # Find the density of each point
            row_density = np.array([self.point_density(p) for p in data[i]])
            # Combine points and normalize
            combined = row_density.sum(axis=0)
            density[i, :] = combined / combined.sum()
        return density


base = BaseMethod(d=D)
# Only running base method on first 2 rows of data because it's slow
density = base.base_user_density(data[:2])
print(density[:2, :5])


class OptimizedMethod:

    def __init__(self, d=50, norm_val_n=50):
        self.d = d
        self.norm_val_n = norm_val_n
        self.bins = np.linspace(0, 1, d + 1).tolist()

        # Generate a standardized set of values to add to each sample to simulate normal distribution
        self.norm_vals = np.array([norm.ppf(x / norm_val_n) * 0.05 for x in range(1, norm_val_n, 1)])

    @timeit
    def optimized_user_density(self, data):

        samples = np.empty((data.shape[0], data.shape[1], self.norm_val_n - 1))
        # transform datapoints to normal distributions around datapoint
        for i in range(self.norm_vals.shape[0]):
            samples[:, :, i] = data + self.norm_vals[i]
        samples = samples.reshape(samples.shape[0], -1)

        #wrap around [0, 1]
        samples = samples % 1

        #loop over samples for density
        density = np.empty((data.shape[0], self.d))
        for i in range(samples.shape[0]):
            hist, bins = np.histogram(samples[i], bins=self.bins)
            density[i, :] = hist / hist.sum()

        return density


om = OptimizedMethod()
#Run optimized method on first 2 rows for apples to apples comparison
density = om.optimized_user_density(data[:2])
#Run optimized method on full data
density = om.optimized_user_density(data)
print(density[:2, :5])

在我的系统上运行，原始方法运行2行数据大约需要8.4秒，而优化方法运行2行数据需要1毫秒，4.7秒完成10000行。我为每种方法打印了前 2 个样本的前五个值。

'base_user_density'  8415.03 ms
[[0.02176227 0.02278653 0.02422535 0.02597123 0.02745976]
 [0.0175103  0.01638513 0.01524853 0.01432158 0.01391156]]
'optimized_user_density'  1.09 ms
'optimized_user_density'  4755.49 ms
[[0.02142857 0.02244898 0.02530612 0.02612245 0.0277551 ]
 [0.01673469 0.01653061 0.01510204 0.01428571 0.01326531]]

Answer 3

通过使用 FFT 并以 numpy 友好格式创建 data，对于最大情况（N=10000，AVG[P]=100，D=50），您可以获得低于 50 毫秒。否则会接近 300 毫秒。

这个想法是将一个以 0 为中心的单一正态分布与一系​​列狄拉克增量进行卷积。

见下图：

使用循环卷积解决了两个问题。

• 自然地处理边缘的包裹
• 可以使用 FFT 和Convolution Theorem 高效计算

第一个必须创建要复制的分发。函数mk_bell() 创建了一个以 0 为中心的标准差 0.05 正态分布的直方图。 分布围绕 1。可以在此处使用 任意 分布。计算出的分布谱用于快速卷积。

接下来创建一个类似梳子的函数。峰值放置在与用户密度峰值相对应的索引处。例如

peaks_location = [0.1, 0.3, 0.7]
D = 10

映射到

peak_index = (D * peak_location).astype(int) = [1, 3, 7]
dist = [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0] # ones at [1, 3, 7]

您可以借助 np.bincount() 函数计算每个峰值位置的 bin 索引，从而快速创建 Diract Delta 的组合。 为了加快速度，可以并行计算用户峰值的梳状函数。

数组dist 是形状为NxD 的二维数组。它可以线性化为形状为(N*D) 的一维数组。在位置[user_id, peak_index] 上的此更改元素之后，将可以从索引user_id*D + peak_index 访问。 使用 numpy 友好的输入格式（如下所述），此操作很容易矢量化。

卷积定理说两个信号的卷积谱等于每个信号的谱的乘积。

频谱是使用numpy.fft.rfft 计算的，这是专用于仅实数信号（无虚部）的快速傅里叶变换的变体。

Numpy 允许使用一个命令计算较大矩阵的每一行的 FFT。

接下来，卷积的频谱通过简单的乘法和广播的使用来计算。

接下来，通过在numpy.fft.irfft 中实现的傅里叶逆变换将频谱计算回“时间”域。

要使用 numpy 的全速，应避免可变大小的数据结构并保持固定大小的数组。我建议将输入数据表示为三个数组。

• uids用户标识，整数0..N-1
• peaks，峰所在位置
• mass，peek 的质量，目前是 1/numer-of-peaks-for-user

这种数据表示允许快速矢量化处理。 例如：

user_data = [[0.1, 0.3], [0.5]]

映射到：

uids = [0, 0, 1] # 2 points for user_data[0], one from user_data[1]
peaks = [0.1, 0.3, 0.5] # serialized user_data
mass = [0.5, 0.5, 1] # scaling factors for each peak, 0.5 means 2 peaks for user 0

代码：

import numpy as np
import matplotlib.pyplot as plt
import time

def mk_bell(D, SIGMA):
    # computes normal distribution wrapped and centered at zero
    x = np.linspace(0, 1, D, endpoint=False);
    x = (x + 0.5) % 1 - 0.5
    bell = np.exp(-0.5*np.square(x / SIGMA))
    return bell / bell.sum()

def user_densities_by_fft(uids, peaks, mass, D, N=None):
    bell = mk_bell(D, 0.05).astype('f4')
    sbell = np.fft.rfft(bell)
    if N is None:
        N = uids.max() + 1
    # ensure that peaks are in [0..1) internal
    peaks = peaks - np.floor(peaks)
    # convert peak location from 0-1 to the indices
    pidx = (D * (peaks + uids)).astype('i4')
    dist = np.bincount(pidx, mass, N * D).reshape(N, D)
    # process all users at once with Convolution Theorem
    sdist = np.fft.rfft(dist)
    sdist *= sbell
    res = np.fft.irfft(sdist)

    return res

def generate_data(N, Pmean):
    # generateor for large data
    data = []
    for n in range(N):
        # select P uniformly from 1..2*Pmean
        P = np.random.randint(2 * Pmean) + 1
        # select peak locations
        chunk = np.random.uniform(size=P)
        data.append(chunk.tolist())
    return data

def make_data_numpy_friendly(data):
    uids = []
    chunks = []
    mass = []
    for uid, peaks in enumerate(data):
        uids.append(np.full(len(peaks), uid))
        mass.append(np.full(len(peaks), 1 / len(peaks)))
        chunks.append(peaks)
    return np.hstack(uids), np.hstack(chunks), np.hstack(mass)




D = 50

# demo for simple multi-distribution
data, N = [[0, .5], [.7, .7, .7, .9], [0.05, 0.3, 0.5, 0.5]], None
uids, peaks, mass = make_data_numpy_friendly(data)
dist = user_densities_by_fft(uids, peaks, mass, D, N)
plt.plot(dist.T)
plt.show()

# the actual measurement
N = 10000
P = 100
data = generate_data(N, P)

tic = time.time()
uids, peaks, mass = make_data_numpy_friendly(data)
toc = time.time()
print(f"make_data_numpy_friendly: {toc - tic}")

tic = time.time()
dist = user_densities_by_fft(uids, peaks, mass, D, N)
toc = time.time()
print(f"user_densities_by_fft: {toc - tic}")

我的 4 核 Haswell 机器上的结果是：

make_data_numpy_friendly: 0.2733159065246582
user_densities_by_fft: 0.04064297676086426

处理数据需要 40 毫秒。请注意，将数据处理为 numpy 友好格式所花费的时间是实际计算分布的 6 倍。 Python 在循环方面真的很慢。 因此，我强烈建议首先以对 numpy 友好的方式直接生成输入数据。

有一些问题需要解决：

• 精度，可以通过使用更大的D 和下采样来提高
• 可以通过加宽尖峰进一步提高峰定位的准确性。
• 性能，scipy.fft 提供可能更快的 FFT 实现的移动变体