46

假设您有一个浮点数列表，这些浮点数大约是一个公共数量的倍数，例如

2.468、3.700、6.1699

大约是 1.234 的所有倍数。你如何描述这个“近似 gcd”，你将如何计算或估计它？

与我对这个问题的回答密切相关。

4

8 回答 8

25

您可以运行 Euclid 的 gcd 算法，任何小于 0.01（或您选择的少量数字）的值都是伪 0。使用您的数字：

3.700 = 1 * 2.468 + 1.232,
2.468 = 2 * 1.232 + 0.004.

所以前两个数字的伪gcd是1.232。现在你用你的最后一个数字来获取这个 gcd：

6.1699 = 5 * 1.232 + 0.0099.

所以1.232是伪gcd，复数是2,3,5。为了改善这个结果，您可以对数据点进行线性回归：

(2,2.468), (3,3.7), (5,6.1699).

斜率是改进后的伪 gcd。

警告：第一部分是算法在数值上不稳定——如果你从非常脏的数据开始，你就有麻烦了。

于 2009-01-16T14:15:38.167 回答
14

将您的测量值表示为最低值的倍数。因此，您的列表变为 1.00000、1.49919、2.49996。这些值的小数部分将非常接近 1/Nths，因为 N 的某个值取决于您的最低值与基频的接近程度。我建议循环增加 N 直到找到足够精细的匹配。在这种情况下，对于 N=1（即假设 X=2.468 是您的基频），您会发现标准偏差为 0.3333（三个值中的两个与 X * 1 相差 0.5），这是不可接受的高。对于 N=2（即假设 2.468/2 是您的基频），您会发现标准偏差几乎为零（所有三个值都在 X/2 的倍数的 0.001 以内），因此 2.468/2 是您的近似值GCD。

我计划中的主要缺陷是，当最低测量值最准确时，它的效果最好，但情况可能并非如此。这可以通过多次执行整个操作来缓解，每次丢弃测量列表中的最小值，然后使用每次通过的结果列表来确定更精确的结果。另一种改进结果的方法是调整 GCD 以最小化 GCD 的整数倍与测量值之间的标准偏差。

于 2009-01-15T00:02:48.737 回答
14

这让我想起了寻找实数的良好有理数近似的问题。标准技术是连续分数展开：

def rationalizations(x):
    assert 0 <= x
    ix = int(x)
    yield ix, 1
    if x == ix: return
    for numer, denom in rationalizations(1.0/(x-ix)):
        yield denom + ix * numer, numer

我们可以将其直接应用于 Jonathan Leffler 和 Sparr 的方法：

>>> a, b, c = 2.468, 3.700, 6.1699
>>> b/a, c/a
(1.4991896272285252, 2.4999594813614263)
>>> list(itertools.islice(rationalizations(b/a), 3))
[(1, 1), (3, 2), (925, 617)]
>>> list(itertools.islice(rationalizations(c/a), 3))
[(2, 1), (5, 2), (30847, 12339)]

从每个序列中挑选出第一个足够好的近似值。（这里是 3/2 和 5/2。）或者不是直接比较 3.0/2.0 和 1.499189...，您可能会注意到 925/617 使用的整数比 3/2大得多，因此 3/2 是一个很好的停止位置.

您除以哪个数字并不重要。（例如，使用 a/b 和 c/b 可以得到 2/3 和 5/3。）一旦有了整数比率，就可以使用 shsmurfy 的线性回归来细化基本面的隐含估计。每个人都赢了！

于 2009-01-26T07:19:14.447 回答
5

我假设你所有的数字都是整数值的倍数。对于我的其余解释，A 将表示您尝试查找的“根”频率，B 将是您必须开始的数字数组。

您正在尝试做的是表面上类似于线性回归。您正在尝试找到一个线性模型 y=mx+b 来最小化线性模型和一组数据之间的平均距离。在您的情况下，b=0，m 是根频率，y 表示给定值。最大的问题是没有明确给出自变量 X。关于 X，我们唯一知道的是它的所有成员都必须是整数。

您的首要任务是尝试确定这些自变量。目前我能想到的最佳方法是假设给定频率具有几乎连续的索引 ( x_1=x_0+n)。所以B_0/B_1=(x_0)/(x_0+n)给定一个（希望是）小整数n。然后，您可以利用这一事实x_0 = n/(B_1-B_0)，从 n=1 开始，并不断提高它，直到 k-rnd(k) 在某个阈值内。在获得 x_0（初始索引）之后，您可以近似根频率 ( A = B_0/x_0)。然后您可以通过查找来近似其他索引x_n = rnd(B_n/A)。这种方法不是很健壮，如果数据中的错误很大，可能会失败。

如果您想要更好地逼近根频率 A，您可以使用线性回归来最小化线性模型的误差，因为您拥有相应的因变量。最简单的方法是使用最小二乘拟合。Wolfram 的数学世界对这个问题进行了深入的数学处理，但通过谷歌搜索可以找到相当简单的解释。

于 2009-01-15T04:31:39.957 回答
4

有趣的问题...不容易。

我想我会看看样本值的比率：

  • 3.700 / 2.468 = 1.499...
  • 6.1699 / 2.468 = 2.4999...
  • 6.1699 / 3.700 = 1.6675...

然后我会在这些结果中寻找一个简单的整数比率。

  • 1.499 ~= 3/2
  • 2.4999 ~= 5/2
  • 1.6675 ~= 5/3

我没有追过它，但是沿着这条线的某个地方，你认为 1:1000 或其他东西的错误已经足够好，然后你回溯以找到基本的近似 GCD。

于 2009-01-15T00:08:50.970 回答
3

我自己看到和使用的解决方案是选择一些常数，比如 1000，将所有数字乘以这个常数，将它们四舍五入，使用标准算法找到这些整数的 GCD，然后将结果除以所述常数(1000)。常数越大，精度越高。

于 2009-01-30T22:53:08.080 回答
1

我在 MathStackExchange（此处和此处）中发现了这个问题，正在寻找我的答案。

我只设法（尚未）测量给定谐波频率列表（遵循声音/音乐命名法）的基频的吸引力，如果您的选项数量减少并且计算吸引力是可行的，这可能很有用每一个，然后选择最合适的。

我在 MSE 的问题中的 C&P（格式更漂亮）：

  • 为 v 列表 {v_1, v_2, ..., v_n}，从低到高排序
  • mean_sin(v, x) = sum(sin(2*pi*v_i/x), for i in {1, ...,n})/n
  • mean_cos(v, x) = sum(cos(2*pi*v_i/x), for i in {1, ...,n})/n
  • gcd_appeal (v, x) = 1 - sqrt(mean_sin(v, x)^2 + (mean_cos(v, x) - 1)^2)/2，产生区间 [0,1] 中的数字。

目标是找到使吸引力最大化的 x。这是您的示例 [2.468, 3.700, 6.1699] 的 ( gcd_appeal ) 图，您发现最佳 GCD 位于x = 1.2337899957639993 在此处输入图像描述

编辑： 您可能会发现这个 JAVA 代码很方便，可以计算除数相对于股息列表的（模糊）可除性（又名 gcd_appeal）；你可以用它来测试你的哪个候选人是最好的除数。代码看起来很丑，因为我试图优化它的性能。

    //returns the mean divisibility of dividend/divisor as a value in the range [0 and 1]
    // 0 means no divisibility at all
    // 1 means full divisibility
    public double divisibility(double divisor, double... dividends) {
        double n = dividends.length;
        double factor = 2.0 / divisor;
        double sum_x = -n;
        double sum_y = 0.0;
        double[] coord = new double[2];
        for (double v : dividends) {
            coordinates(v * factor, coord);
            sum_x += coord[0];
            sum_y += coord[1];
        }
        double err = 1.0 - Math.sqrt(sum_x * sum_x + sum_y * sum_y) / (2.0 * n);
        //Might happen due to approximation error
        return err >= 0.0 ? err : 0.0;
    }

    private void coordinates(double x, double[] out) {
        //Bhaskara performant approximation to
        //out[0] = Math.cos(Math.PI*x);
        //out[1] = Math.sin(Math.PI*x);
        long cos_int_part = (long) (x + 0.5);
        long sin_int_part = (long) x;
        double rem = x - cos_int_part;
        if (cos_int_part != sin_int_part) {
            double common_s = 4.0 * rem;
            double cos_rem_s = common_s * rem - 1.0;
            double sin_rem_s = cos_rem_s + common_s + 1.0;
            out[0] = (((cos_int_part & 1L) * 8L - 4L) * cos_rem_s) / (cos_rem_s + 5.0);
            out[1] = (((sin_int_part & 1L) * 8L - 4L) * sin_rem_s) / (sin_rem_s + 5.0);
        } else {
            double common_s = 4.0 * rem - 4.0;
            double sin_rem_s = common_s * rem;
            double cos_rem_s = sin_rem_s + common_s + 3.0;
            double common_2 = ((cos_int_part & 1L) * 8L - 4L);
            out[0] = (common_2 * cos_rem_s) / (cos_rem_s + 5.0);
            out[1] = (common_2 * sin_rem_s) / (sin_rem_s + 5.0);
        }
    }
于 2014-09-01T16:19:08.357 回答
1

当您先验地选择 3 个正公差 (e1,e2,e3) 时，这是对 shsmurfy 解决方案的重新制定，
然后问题是搜索最小的正整数 (n1,n2,n3) 并因此搜索最大的根频率 f，使得：

f1 = n1*f +/- e1
f2 = n2*f +/- e2
f3 = n3*f +/- e3

我们假设 0 <= f1 <= f2 <= f3
如果我们固定 n1，那么我们得到以下关系：

f  is in interval I1=[(f1-e1)/n1 , (f1+e1)/n1]
n2 is in interval I2=[n1*(f2-e2)/(f1+e1) , n1*(f2+e2)/(f1-e1)]
n3 is in interval I3=[n1*(f3-e3)/(f1+e1) , n1*(f3+e3)/(f1-e1)]

我们从 n1 = 1 开始，然后递增 n1 直到区间 I2 和 I3 包含一个整数 - 这floor(I2min) different from floor(I2max)与 I3 相同
然后我们在区间 I2 中选择最小整数 n2，在区间 I3 中选择最小整数 n3。

假设浮点误差的正态分布，根频率 f 的最可能估计是最小化的

J = (f1/n1 - f)^2 + (f2/n2 - f)^2 + (f3/n3 - f)^2

那是

f = (f1/n1 + f2/n2 + f3/n3)/3

如果在区间 I2,I3 中有多个整数 n2,n3，我们也可以选择最小化残差的对

min(J)*3/2=(f1/n1)^2+(f2/n2)^2+(f3/n3)^2-(f1/n1)*(f2/n2)-(f1/n1)*(f3/n3)-(f2/n2)*(f3/n3)

另一种变体可能是继续迭代并尝试最小化另一个标准，如 min(J(n1))*n1，直到 f 低于某个频率（n1 达到上限）......

于 2012-11-15T00:27:36.513 回答