【问题标题】:How can I calculate the point between two overlapping linear datasets?如何计算两个重叠线性数据集之间的点?
【发布时间】:2014-09-29 23:09:33
【问题描述】:

我有两组数据有些重叠(见下图)。我需要找到这些集合之间的点,在那里人们会猜测未知数据点属于特定类别。

如果我有一个新的数据点(比如5000),并且必须赌 $$$ 它属于 A 组还是 B 组,我该如何计算该点让我的赌注最有把握?

请参阅下面的示例数据集和随附图表,以及这些组之间的近似点(通过眼睛计算)。

GROUP A
[385,515,975,1136,2394,2436,4051,4399,4484,4768,4768,4849,4856,4954,5020,5020,5020,5020,5020,5020,5020,5020,5020,5052,5163,5200,5271,5421,5421,5442,5746,5765,5903,5992,5992,6046,6122,6205,6208,6239,6310,6360,6416,6512,6536,6543,6581,6609,6696,6699,6752,6796,6806,6855,6859,6886,6906,6911,6923,6953,7016,7072,7086,7089,7110,7232,7278,7293,7304,7309,7348,7367,7378,7380,7419,7453,7454,7492,7506,7549,7563,7721,7723,7731,7745,7750,7751,7783,7791,7813,7813,7814,7818,7833,7863,7875,7886,7887,7902,7907,7935,7942,7942,7948,7973,7995,8002,8013,8013,8015,8024,8025,8030,8038,8041,8050,8056,8060,8064,8071,8081,8082,8085,8093,8124,8139,8142,8167,8179,8204,8214,8223,8225,8247,8248,8253,8258,8264,8265,8265,8269,8277,8278,8289,8300,8312,8314,8323,8328,8334,8363,8369,8390,8397,8399,8399,8401,8436,8442,8456,8457,8471,8474,8483,8503,8511,8516,8533,8560,8571,8575,8583,8592,8593,8626,8635,8635,8644,8659,8685,8695,8695,8702,8714,8715,8717,8729,8732,8740,8743,8750,8756,8772,8772,8778,8797,8828,8840,8840,8843,8856,8865,8874,8876,8878,8885,8887,8893,8896,8905,8910,8955,8970,8971,8991,8995,9014,9016,9042,9043,9063,9069,9104,9106,9107,9116,9131,9157,9227,9359,9471]

GROUP B
[12,16,29,32,33,35,39,42,44,44,44,45,45,45,45,45,45,45,45,45,47,51,51,51,57,57,60,61,61,62,71,75,75,75,75,75,75,76,76,76,76,76,76,79,84,84,85,89,93,93,95,96,97,98,100,100,100,100,100,102,102,103,105,108,109,109,109,109,109,109,109,109,109,109,109,109,110,110,112,113,114,114,116,116,118,119,120,121,122,124,125,128,129,130,131,132,133,133,137,138,144,144,146,146,146,148,149,149,150,150,150,151,153,155,157,159,164,164,164,167,169,170,171,171,171,171,173,174,175,176,176,177,178,179,180,181,181,183,184,185,187,191,193,199,203,203,205,205,206,212,213,214,214,219,224,224,224,225,225,226,227,227,228,231,234,234,235,237,240,244,245,245,246,246,246,248,249,250,250,251,255,255,257,264,264,267,270,271,271,281,282,286,286,291,291,292,292,294,295,299,301,302,304,304,304,304,304,306,308,314,318,329,340,344,345,356,359,363,368,368,371,375,379,386,389,390,392,394,408,418,438,440,456,456,458,460,461,467,491,503,505,508,524,557,558,568,591,609,622,656,665,668,687,705,728,817,839,965,1013,1093,1126,1512,1935,2159,2384,2424,2426,2484,2738,2746,2751,3006,3184,3184,3184,3184,3184,4023,5842,5842,6502,7443,7781,8132,8237,8501]

数组统计:

                      Group A    Group B
Total Numbers             231        286
Mean                  7534.71     575.56
Standard Deviation    1595.04    1316.03

【问题讨论】:

  • 你知道A组/B组的分布吗?

标签: arrays algorithm sorting machine-learning classification


【解决方案1】:

在合理的假设下,一个好的判别式是唯一的数据值,它导致 B 的概率密度在分割点左侧的面积等于 A 在右侧的面积(反之亦然,给出相同的点) .

一个简单的方法是计算两个经验累积分布函数 (CDF) 作为shown here 并搜索它们以提供分割点。这是两个 CDF 之和为 1 的点。

简单地说,构建经验 CDF 只是对每个数据集进行排序并将数据用作 x 轴值。从左到右绘制曲线,从 y=0 开始,在每个 x 值处向上走 1/n 步。这样的曲线从 x 1 的 0 逐渐上升到 x >= data[n] 的 y = CDF(x) = 1。有一种稍微复杂的方法可以给出连续的逐步线性曲线而不是阶梯,在某些假设下,它是对真实 CDF 的更好估计。在

请注意,上面的讨论只是为了提供直觉。 CDF 由已排序的数据数组完美表示。不需要新的数据结构;即 x[i], i=1,2,...,n 是曲线到达 y = i/n 处的 x 值。

使用两个 CDF,R(x) 和 B(x),根据您的图表,您想要找到唯一点 x,使得 |1 - R(x) - B(x)|被最小化(使用分段线性 CDF,您将始终能够将其设为零)。这可以通过二分搜索很容易地完成。

这种方法的一个好处是您可以通过将两个 CDF 维护在有序集合(平衡二叉搜索树)中来使其动态化。随着点的增加,新的分割点很容易找到。

有序集需要“订单统计”。 Here is a reference。我的意思是您需要能够查询排序集以检索 CDF 中任何存储的 x 值的序数。这可以通过跳过列表和树来完成。

我编写了这个算法的一个变体。它使用分段 CDF 近似,但也允许在重复数据点处进行“垂直步骤”。这使算法有些复杂,但还不错。然后我使用二分法(而不是组合二分搜索)来找到分割点。正常的二分算法需要修改以适应 CDF 中的垂直“步骤”。我认为这一切都是正确的,但经过了轻微的测试。

处理的一个极端情况是数据集具有不相交的范围。这将在较低的顶部和较高的底部之间找到 a 点,这是一个完全有效的鉴别器。但是你可能想做一些更有趣的事情,比如返回某种加权平均值。

请注意,如果您对数据可以达到的实际最小值和最大值有一个很好的概念,并且它们不会出现在数据中,则应考虑添加它们,以免 CDF不经意间偏了。

在您的示例数据上,代码生成 4184.76,看起来非常接近您在图表中选择的值(略低于最小和最大数据的一半)。

请注意,我没有对数据进行排序,因为它已经是。排序绝对是必要的。

public class SplitData {

    // Return: i such that a[i] <= x < a[i+1] if i,i+1 in range
    // else -1 if x < a[0]
    // else a.length if x >= a[a.length - 1]
    static int hi_bracket(double[] a, double x) {
        if (x < a[0]) return -1;
        if (x >= a[a.length - 1]) return a.length;
        int lo = 0, hi = a.length - 1;
        while (lo + 1 < hi) {
            int mid = (lo + hi) / 2;
            if (x < a[mid])
                hi = mid;
            else 
                lo = mid;
        }
        return lo;
    }

    // Return: i such that a[i-1] < x <= a[i] if i-1,i in range
    // else -1 if x <= a[0]
    // else a.length if x > a[a.length - 1]
    static int lo_bracket(double[] a, double x) {
        if (x <= a[0]) return -1;
        if (x > a[a.length - 1]) return a.length;
        int lo = 0, hi = a.length - 1;
        while (lo + 1 < hi) {
            int mid = (lo + hi) / 2;
            if (x <= a[mid])
                hi = mid;
            else
                lo = mid;
        }
        return hi;
    }

    // Interpolate the CDF value for the data a at value x.  Returns a range.
    static void interpolate_cdf(double[] a, double x, double[] rtn) {
        int lo_i1 = lo_bracket(a, x);
        if (lo_i1 == -1) {
            rtn[0] = rtn[1] = 0;
            return;
        }
        int hi_i0 = hi_bracket(a, x);
        if (hi_i0 == a.length) {
            rtn[0] = rtn[1] = 1;
            return;
        }
        if (hi_i0 + 1 == lo_i1) {  // normal interpolation
            rtn[0] = rtn[1] 
                = (hi_i0 + (x - a[hi_i0]) / (a[lo_i1] - a[hi_i0])) 
                    / (a.length - 1);
            return;
        }
        // we're on a joint or step; return range answer
        rtn[0] = (double)lo_i1 / (a.length - 1);  
        rtn[1] = (double)hi_i0 / (a.length - 1);
        assert rtn[0] <= rtn[1];
    }

    // Find the data value where the two given data set's empirical CDFs
    // sum to 1. This is a good discrimination value for new data.
    // This deals with the case where there's a step in either or both CDFs.
    static double find_bisector(double[] a, double[] b) {
        assert a.length > 0;
        assert b.length > 0;
        double lo = Math.min(a[0], b[0]);
        double hi = Math.max(a[a.length - 1], b[b.length - 1]);
        double eps = (hi - lo) * 1e-7;
        double[] a_rtn = new double[2], b_rtn = new double[2];
        while (hi - lo > eps) {
            double mid = 0.5 * (lo + hi);
            interpolate_cdf(a, mid, a_rtn);
            interpolate_cdf(b, mid, b_rtn);
            if (1 < a_rtn[0] + b_rtn[0])
                hi = mid;
            else if (a_rtn[1] + b_rtn[1] < 1)
                lo = mid;
            else
                return mid;  // 1 is included in the interpolated range
        }
        return 0.5 * (lo + hi);
    }

    public static void main(String[] args) {
        double split = find_bisector(a, b);
        System.err.println("Split at x = " + split);
    }

    static final double[] a = {
        385, 515, 975, 1136, 2394, 2436, 4051, 4399, 4484, 4768, 4768, 4849,
        4856, 4954, 5020, 5020, 5020, 5020, 5020, 5020, 5020, 5020, 5020, 5052,
        5163, 5200, 5271, 5421, 5421, 5442, 5746, 5765, 5903, 5992, 5992, 6046,
        6122, 6205, 6208, 6239, 6310, 6360, 6416, 6512, 6536, 6543, 6581, 6609,
        6696, 6699, 6752, 6796, 6806, 6855, 6859, 6886, 6906, 6911, 6923, 6953,
        7016, 7072, 7086, 7089, 7110, 7232, 7278, 7293, 7304, 7309, 7348, 7367,
        7378, 7380, 7419, 7453, 7454, 7492, 7506, 7549, 7563, 7721, 7723, 7731,
        7745, 7750, 7751, 7783, 7791, 7813, 7813, 7814, 7818, 7833, 7863, 7875,
        7886, 7887, 7902, 7907, 7935, 7942, 7942, 7948, 7973, 7995, 8002, 8013,
        8013, 8015, 8024, 8025, 8030, 8038, 8041, 8050, 8056, 8060, 8064, 8071,
        8081, 8082, 8085, 8093, 8124, 8139, 8142, 8167, 8179, 8204, 8214, 8223,
        8225, 8247, 8248, 8253, 8258, 8264, 8265, 8265, 8269, 8277, 8278, 8289,
        8300, 8312, 8314, 8323, 8328, 8334, 8363, 8369, 8390, 8397, 8399, 8399,
        8401, 8436, 8442, 8456, 8457, 8471, 8474, 8483, 8503, 8511, 8516, 8533,
        8560, 8571, 8575, 8583, 8592, 8593, 8626, 8635, 8635, 8644, 8659, 8685,
        8695, 8695, 8702, 8714, 8715, 8717, 8729, 8732, 8740, 8743, 8750, 8756,
        8772, 8772, 8778, 8797, 8828, 8840, 8840, 8843, 8856, 8865, 8874, 8876,
        8878, 8885, 8887, 8893, 8896, 8905, 8910, 8955, 8970, 8971, 8991, 8995,
        9014, 9016, 9042, 9043, 9063, 9069, 9104, 9106, 9107, 9116, 9131, 9157,
        9227, 9359, 9471
    };
    static final double[] b = {
        12, 16, 29, 32, 33, 35, 39, 42, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45,
        45, 45, 47, 51, 51, 51, 57, 57, 60, 61, 61, 62, 71, 75, 75, 75, 75, 75,
        75, 76, 76, 76, 76, 76, 76, 79, 84, 84, 85, 89, 93, 93, 95, 96, 97, 98,
        100, 100, 100, 100, 100, 102, 102, 103, 105, 108, 109, 109, 109, 109,
        109, 109, 109, 109, 109, 109, 109, 109, 110, 110, 112, 113, 114, 114,
        116, 116, 118, 119, 120, 121, 122, 124, 125, 128, 129, 130, 131, 132,
        133, 133, 137, 138, 144, 144, 146, 146, 146, 148, 149, 149, 150, 150,
        150, 151, 153, 155, 157, 159, 164, 164, 164, 167, 169, 170, 171, 171,
        171, 171, 173, 174, 175, 176, 176, 177, 178, 179, 180, 181, 181, 183,
        184, 185, 187, 191, 193, 199, 203, 203, 205, 205, 206, 212, 213, 214,
        214, 219, 224, 224, 224, 225, 225, 226, 227, 227, 228, 231, 234, 234,
        235, 237, 240, 244, 245, 245, 246, 246, 246, 248, 249, 250, 250, 251,
        255, 255, 257, 264, 264, 267, 270, 271, 271, 281, 282, 286, 286, 291,
        291, 292, 292, 294, 295, 299, 301, 302, 304, 304, 304, 304, 304, 306,
        308, 314, 318, 329, 340, 344, 345, 356, 359, 363, 368, 368, 371, 375,
        379, 386, 389, 390, 392, 394, 408, 418, 438, 440, 456, 456, 458, 460,
        461, 467, 491, 503, 505, 508, 524, 557, 558, 568, 591, 609, 622, 656,
        665, 668, 687, 705, 728, 817, 839, 965, 1013, 1093, 1126, 1512, 1935,
        2159, 2384, 2424, 2426, 2484, 2738, 2746, 2751, 3006, 3184, 3184, 3184,
        3184, 3184, 4023, 5842, 5842, 6502, 7443, 7781, 8132, 8237, 8501
    };
}

【讨论】:

  • 引人入胜的方法。是的!这很有趣。本周初我将专注于解决这个问题。谢谢。
  • @Ryan 好的。语言?我可以轻松地做 Java、Ruby、C/C++。如果你真的需要 Python(不是我的初选之一)。
  • 我最终会移植到 PHP,所以只要你的方法不依赖于独特的库,我应该也可以。即使是伪代码也可以。
  • @Ryan 对不起,还有一个问题。数据总是整数吗?需要知道 CDF 应该是离散的还是连续的。
  • 总是数字,但不一定是整数。但如有必要,我可以转换为舍入值并退出它。
【解决方案2】:

我只想指出另一种使用密度估计的方法。

根据您的数据,使用kernel density estimation 很容易拟合平滑的pdf。下面的python代码展示了如何使用scipy中的kde模块。

from scipy.stats.kde import gaussian_kde
from numpy import linspace
import matplotlib.pyplot as plt
data1 = [385,515,975,1136,2394,2436,4051,4399,4484,4768,4768,4849,4856,4954,5020,5020,5020,5020,5020,5020,5020,5020,5020,5052,5163,5200,5271,5421,5421,5442,5746,5765,5903,5992,5992,6046,6122,6205,6208,6239,6310,6360,6416,6512,6536,6543,6581,6609,6696,6699,6752,6796,6806,6855,6859,6886,6906,6911,6923,6953,7016,7072,7086,7089,7110,7232,7278,7293,7304,7309,7348,7367,7378,7380,7419,7453,7454,7492,7506,7549,7563,7721,7723,7731,7745,7750,7751,7783,7791,7813,7813,7814,7818,7833,7863,7875,7886,7887,7902,7907,7935,7942,7942,7948,7973,7995,8002,8013,8013,8015,8024,8025,8030,8038,8041,8050,8056,8060,8064,8071,8081,8082,8085,8093,8124,8139,8142,8167,8179,8204,8214,8223,8225,8247,8248,8253,8258,8264,8265,8265,8269,8277,8278,8289,8300,8312,8314,8323,8328,8334,8363,8369,8390,8397,8399,8399,8401,8436,8442,8456,8457,8471,8474,8483,8503,8511,8516,8533,8560,8571,8575,8583,8592,8593,8626,8635,8635,8644,8659,8685,8695,8695,8702,8714,8715,8717,8729,8732,8740,8743,8750,8756,8772,8772,8778,8797,8828,8840,8840,8843,8856,8865,8874,8876,8878,8885,8887,8893,8896,8905,8910,8955,8970,8971,8991,8995,9014,9016,9042,9043,9063,9069,9104,9106,9107,9116,9131,9157,9227,9359,9471]
data2 = [12,16,29,32,33,35,39,42,44,44,44,45,45,45,45,45,45,45,45,45,47,51,51,51,57,57,60,61,61,62,71,75,75,75,75,75,75,76,76,76,76,76,76,79,84,84,85,89,93,93,95,96,97,98,100,100,100,100,100,102,102,103,105,108,109,109,109,109,109,109,109,109,109,109,109,109,110,110,112,113,114,114,116,116,118,119,120,121,122,124,125,128,129,130,131,132,133,133,137,138,144,144,146,146,146,148,149,149,150,150,150,151,153,155,157,159,164,164,164,167,169,170,171,171,171,171,173,174,175,176,176,177,178,179,180,181,181,183,184,185,187,191,193,199,203,203,205,205,206,212,213,214,214,219,224,224,224,225,225,226,227,227,228,231,234,234,235,237,240,244,245,245,246,246,246,248,249,250,250,251,255,255,257,264,264,267,270,271,271,281,282,286,286,291,291,292,292,294,295,299,301,302,304,304,304,304,304,306,308,314,318,329,340,344,345,356,359,363,368,368,371,375,379,386,389,390,392,394,408,418,438,440,456,456,458,460,461,467,491,503,505,508,524,557,558,568,591,609,622,656,665,668,687,705,728,817,839,965,1013,1093,1126,1512,1935,2159,2384,2424,2426,2484,2738,2746,2751,3006,3184,3184,3184,3184,3184,4023,5842,5842,6502,7443,7781,8132,8237,8501]

pdf1 = gaussian_kde(data1)
pdf2 = gaussian_kde(data2)

x = linspace(0, 9500, 1000)
plt.plot(x, pdf1(x),'r')
plt.plot(x, pdf2(x),'g')
plt.legend(['data1 pdf', 'data2 pdf'])

plt.show()

图中,绿色为第二个数据集的pdf;红色是第一个数据集的 pdf。显然,决策边界是通过绿色与红色相交点的垂直线。

要以数字方式找到边界,我们可以执行如下操作(假设只有一个交叉点,否则没有意义):

min_diff = 10000
min_diff_x = -1
for x in linspace(3600, 4000, 400):
    diff = abs(pdf1(x) - pdf2(x))
    if diff < min_diff:
        min_diff = diff
        min_diff_x = x
print min_diff, min_diff_x

我们发现边界大约位于 3762。

如果两个 pdf 有多个交集,为了预测数据点 x 属于哪个类,我们计算 pdf1(x)pdf2(x),最大的类是最小化贝叶斯风险的类。有关贝叶斯风险和预测错误概率评估的更多详细信息,请参阅here

下图是一个实际包含三个pdf的例子,在任意查询点x,我们应该分别询问这三个pdf,并选择最大值为pdf(x)的那个作为预测类。

【讨论】:

  • 我喜欢使用密度的想法,但我不遵循“清楚”。如果密度更频繁地相交怎么办?也许您可以使用累积分布?
  • @Teepeemm 我的看法是你不想要密度的交集。您想要右侧区域与左侧区域匹配的点,这是唯一的。最好在 CDF 中找到这一点。
  • 当然,如果两个 pdf 相交更多,给定一个新数据点 x,我们只需查看 pdf1(x)pdf2(x) 并将最大值作为我们的预测。而已。我试图找出交叉点只是因为给定的数据集分离良好并且可以提供清晰的决策边界。
  • 我在这里尝试了许多答案,并且能够通过这种方法始终如一地达到预期的结果。谢谢。
【解决方案3】:

您正在描述一个一维statistical classification 问题,您正在寻找“决策边界”。您有很多选择:

  • 逻辑回归
  • 最近邻分类器
  • 支持向量机
  • 多层感知器
  • ...

但由于问题很简单(一维,两个分离良好的类)并且决策边界是一个相当空的区域,我怀疑没有任何繁重的统计方法会显着优于简单的基于眼睛的猜测。

【讨论】:

  • 是的,但重点是避免需要基于眼睛的猜测。
【解决方案4】:

这可以被视为具有单个连续预测器的二元分类问题。您可以将此视为拟合一个简单的决策树,找到一个阈值 t,以便您在值 >= t 时预测 A 组。

为此,您选择最小化结果拆分的熵的 t。假设您对某些 t 有以下计数:

| | <t | >= t | | Group A | X | Y | | Group B | Z | W |

= 分裂的熵是 -(Y/(Y+W))*log(Y/(Y+W)) - (W/(Y+W))*log(W/(Y+W)) .这看起来比实际更混乱;它只是拆分中每个组的比例 p 的 -p*log(p) 的总和。

您取两者的加权平均值,按拆分的整体大小加权。所以第一项由 (X+Z)/(X+Y+Z+W) 加权,另一项由 (Y+W)/(X+Y+Z+W) 加权。​​

【讨论】:

【解决方案5】:

您可以针对每个集合计算新点的Mahalanobis distance。与新点距离最短的集合最有可能匹配。

马氏距离是对点 P 和分布 D 之间距离的度量,由 PC Mahalanobis 在 1936 年引入。1 它是测量距离 P 有多少标准差的思想的多维概括从 D 的平均值。如果 P 处于 D 的平均值,这个距离为零,并且随着 P 远离平均值而增长

由于您的空间是一维的,因此计算应简化为:

  1. 计算每个分布的标准差
  2. 计算每个分布的平均值
  3. 对于每个分布,计算该点与分布均值相差多少标准差。

【讨论】:

  • 谢谢埃里克。那么向量计划将是远离每个分布均值的最大标准差吗?如果 A 组的分数不成比例,这种情况会改变吗?例如,如果 A 组的点数是 B 组的 100 倍,这仍然有效吗?
  • 假设您在 B 中仍有足够的点数可以很好地估计 B 的均值和标准差,它仍然可以工作。
  • 如果您的样本中有统计上有效的点数,它将起作用。这是一个概率估计。样本量越大,确定性就越高。
  • @Ryan 分界线将是在任一组中的概率相同的位置,或者 A 组平均值的标准差数等于 B 组平均值的点。我为您的数据集获得了大约 3721.65
  • @Ryan 我简单地使用了Eric提出的方法,在步骤3中,设置(|meanA - x| / sigmaA) = (|meanB - x| / sigmaB)并求解x
猜你喜欢
  • 2021-03-01
  • 1970-01-01
  • 1970-01-01
  • 2020-05-08
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2016-11-11
  • 2015-01-25
相关资源
最近更新 更多