65.9K
CodeProject 正在变化。 阅读更多。
Home

通过组合数据点拼接实现呼叫分配参数

Fedor Naumenko
starIconstarIconstarIconstarIconstarIcon

5.00/5 (3投票s)

2021年3月29日

CPOL

6分钟阅读

viewsIcon

7412

downloadIcon

110

从分布式物理值集合中查找分布类型和参数

引言

想象一下，我们有一组物理值，它们根据对数正态分布、正态分布或伽马分布进行分布。
我们需要确定分布的类型并调用（找出）其参数。

背景

首先，我们构建值的频率分布。其次，通过分布的图形视图，我们确定其参数。最后，我们将原始分布与计算出的每个分布进行比较，并选择最佳匹配。

仔细观察这些点

  1. 频率分布（概率分布）以典型方式使用排序的关联容器（例如，C++中的std::map）进行构建。

  2. 也许有不止一种方法可以通过外观确定分布参数。我只想到了一种——通过“特征点”的坐标。

    所有分布都属于两参数族。对于所有这些分布，除了主要的概率密度函数之外，这两个参数还通过众数（出现频率最高的值）方程相关联。因此，两个点足以求解两个方程组。其中一个点很明显——这是曲线的顶点，对应于众数。第二个点应属于曲线的上升（左）或下降（右）分支，位于对分布参数最敏感的位置。那里有半高点。其中，最好选择右半高点，因为对于对数正态分布和伽马分布，下降分支通常更平坦，这使得在横坐标上的测量误差更小。

    我们将众数的x坐标表示为M，将半高点（右侧）的x坐标表示为X。
    因此，对于正态分布，我们有

     和  µ = M,

    由此推导出σ方程

    对于对数正态分布，我们有

     和  

    最终得出

    对于伽马分布，我们有

     和  

    最终得出

  3. 两条曲线的比较是皮尔逊方法普遍接受的问题。

概念上看起来很简单。

不幸的是，世界是一个肮脏、混乱、复杂的地方。由于一些物理和技术原因，实际分布经常与理想分布相去甚远。此外，它们经常带有噪声，有时还被调制。更重要的是，它们可能会被截断或在左侧失真（参见DNA片段分布和DNA片段读取分布；所有提供的数据都是真实的）。嗯，通过选择下降分支上的第二个点可以消除左侧的失真和截断。但对于所有其他情况，找到特征点的“真实”坐标是一项非平凡的任务。

问题在于构建一个“平滑”的虚拟曲线，该曲线根据特定标准最能拟合数据，即虚拟点与其真实原型之间的等距（“平滑”意味着沿整个曲线的虚拟点之间的总距离最小）。通常有三种方法可以解决这个问题：非线性回归、贝塞尔曲线和移动平均。现在我们来考虑简化的特殊任务：（a）每个真实点都有一个唯一的横坐标，（b）在包含特征点的关键区域，真实点沿x轴均匀分布，（c）每个虚拟点对应一个真实点（具有相同的横坐标，但经过修正的纵坐标）。在这些条件下，回归和贝塞尔分析都显得过于复杂，因此资源消耗巨大。此外，在调制分布的情况下，贝塞尔曲线无能为力。最佳解决方案似乎是使用基于滚动平均的样条插值方法。

简单移动平均能够完美地完成任务，同时在移动窗口内的真实纵坐标与虚拟纵坐标的总偏差不超过相邻虚拟点纵坐标之差。然而，在存在强局部偏差（“峰值”）的情况下，该方法将无能为力。能够很好地消除峰值的是移动中位数，尽管它实际上并没有对数据进行平均。

解决方案是一致地结合这两种方法。在每次窗口移动时，将预先计算的中位数输入到平均器中。图显示了每个拼接以及组合拼接的结果，分别用于高度噪声数据（a）和光滑但深度调制的（b）。
有趣的是，在“前瞻”窗口（从平均窗口结束的点开始）中计算中位数的变体与在同一窗口中进行群集产生相同的结果。但是，它需要对序列的末尾进行额外的控制（假设Push()方法的输入是迭代器而不是值）。因此，它被放弃了。

现在需要解决移动（滑动）窗口长度的问题。窗口越小，平均峰值纵坐标与“真实”值之间的偏差就越小，从而可以更准确地确定半高点（第二个特征点的纵坐标）。另一方面，对于不光滑或调制的数据，小窗口长度可能导致“虚假”的局部极值。这会导致众数位置不准确，进而导致分布近似不准确（参见图）。该问题分两个阶段解决：首先，根据数据的性质确定最大可靠窗口长度，然后将其循环减小，并在每一步计算皮尔逊相关系数（PCC）。当达到最小窗口长度或连续三步PCC小于最大值时，周期结束。实验表明，这已足够——在PCC连续下降三步后，随着窗口长度的减小，PCC再也不会上升。最大PCC为我们提供了应用分布与真实分布最佳近似的参数。通过应用所有三个定律，我们可以最终回答给定真实分布最可能的“真实”性质——正态、对数正态还是伽马。

拼接器的实现

移动平均和移动中位数拼接器被实现为简单的类，它们继承自移动窗口并具有一个`Push()`方法。窗口长度始终为奇数。构造函数接受一半的长度，因此对于最小半长1，最小窗口长度为3。

// Method declarations and definitions are merged
// solely for the purpose of shortening the presented code

// Moving window (sliding subset)
class MW : protected vector<ULONG>
{
protected:
    // Constructor by moving window half-length ("base")
    MW(size_t base) { insert(begin(), Size(base), 0); }

    // Adds last value and pops the first one (QUEUE functionality)
    void PushVal(ULONG x) {
        move(begin() + 1, end(), begin());
        *--end() = x;
    }

public:
    // Returns length of moving window
    //   @base: moving window half-length
    static size_t Size(size_t base) { return 2 * base + 1; }
};

// Simple Moving Average splicer
class MA : public MW
{
    ULONG _sum = 0;                    // sum of adding values
public:
    // Constructor by moving window half-length ("base")
    inline MA(size_t base) : MW(base) {}

    // Adds value and returns average
    float Push(ULONG x) {
        _sum += x - *begin();
        PushVal(x);
        return float(_sum) / size();
    }
};

// Simple Moving Median splicer
class MM : public MW
{
    vector<ULONG> _ss;                 // sorted sliding subset
public:
    // Constructor by moving window half-length
    MM(size_t base) : MW(base) { _ss.insert(_ss.begin(), size(), 0); }

    // Adds value and returns median
    ULONG Push(ULONG x) {
        PushVal(x);
        copy(begin(), end(), _ss.begin());
        sort(_ss.begin(), _ss.end());
        return _ss[size() >> 1];        // mid-size
    }
};

在遍历真实点的两个运算符中，可以提取平均（虚拟）点的坐标

// Defines key points
//   @base: moving window half-length
//   @summit: returned X,Y coordinates of spliced (smoothed) summit
//   return: key points: X-coord of highest point, X-coord of right middle hight point
fpair LenFreq::GetKeyPoints(size_t base, point& summit) const
{
    point p0(*begin()), p(0, 0);         // previous, current spliced point
    MA ma(base);
    MM mm(base);

    summit.second = 0;
    for (const value_type& f : *this) {  // loop on real points
        p.first = f.first – 2 * base;               // X: minus MA & MM base back shift
        p.second = ma.Push( mm.Push(f.second) );    // Y: spliced
        if (p.second >= summit.second)    summit = p;
    }

    return fpair(
        float(summit.first),             // summit X
        // final point with half height; proportional X
        p0.first + p0.second / (p.second + p0.second)    
    );
}

此处应用了一个稍简化的测试版本。完整实现的callDist实用程序可以在GitHub上找到。

效率

在2.5 GHz笔记本电脑上测量了调用对数正态参数的运行时间。
下表显示了整个算法的最佳和最差运行时间，平均重复1000次。

数据（图示见片段分布） 真实点数 达到最佳窗口长度的迭代次数 最佳窗口长度 时间，微秒
19 2985 2 3（最小值） 46
11 8730 2 3（最小值） 80
17* 5600 6 39 2500

* 本文插图中使用的“噪声”数据

致谢

感谢Oleg Tolmachev（英国格林威治大学）协助改进本文的语言。

历史

  • 2021年3月29日：初始版本
© . All rights reserved.