“近似”最大公约数

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

2.468,3.700,6.1699

这几乎是1.234的整数倍。 你如何描述这个“近似的gcd”,你将如何继续计算或估计它?

与我对这个问题的回答严格相关。

你可以运行Euclid的gcdalgorithm,任何小于0.01(或者你select的一小部分)都是伪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。

警告:这个algorithm的第一部分是algorithm在数值上是不稳定的 – 如果你从非常脏的数据开始,你有麻烦。

将您的测量结果显示为最低点的倍数。 因此你的清单变成了1.00000,1.49919,2.49996。 这些值的小数部分将非常接近1 / N,N的某个值由最低值与基本频率的接近程度决定。 我会build议循环增加N,直到find足够精确的匹配。 在这种情况下,对于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的整数倍与测量值之间的标准偏差最小化。

这让我想起寻找实数有理数逼近的问题。 标准技术是一个持续的分数扩张:

 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的线性回归来优化基本的隐含估计。 每个人都赢了!

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

你试图做的是表面上类似于线性回归 。 您正试图find一个线性模型y = mx + b,它使线性模型与一组数据之间的平均距离最小化。 在你的情况下,b = 0,m是根频率,y代表给定的值。 最大的问题是自variablesX没有明确给出。 我们唯一知道的X是它的所有成员必须是整数。

你的第一个任务是试图确定这些自variables。 我现在所能想到的最好的方法是假定给定的频率具有几乎连续的索引( 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,则可以使用线性回归将线性模型的误差减至最小,现在您有相应的因variables。 最简单的方法是使用最小二乘拟合。 Wolfram的Mathworld对这个问题有一个深入的math处理,但是一些简单的解释可以通过一些Googlesearchfind。

有趣的问题…不容易。

我想我会看样本值的比例:

  • 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或者某个东西的错误已经足够好了,然后你回过头来find基本近似的GCD。

我所看到和使用过的解决scheme是select一些常数,比如1000,用这个常数乘以所有数字,将它们四舍五入到整数,用标准algorithmfind这些整数的GCD,然后用上述常数(1000)。 常数越大,精度越高。

当你先验select3个正面容忍度(e1,e2,e3)时,这是一个shsmurfy解决scheme的改革。
那么问题是search最小的正整数(n1,n2,n3),从而search最大的根频率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的间隔包含一个整数 – 即与I2 floor(I2min) different from floor(I2max)
然后我们select区间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,我们也可以select使剩余最小化的一对

 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达到上限)…

我发现这个问题在MathStackExchange( 这里和这里 )寻找我的答案。

我只pipe理(还)测量一个基本频率的吸引力 ,给定一个谐波频率列表(遵循声音/音乐命名法),如果您的选项数量减less并且可以计算出上诉每一个然后select最适合。

C&P从我的MSE问题(格式更漂亮):

  • 成为列表{v_1,v_2,…,v_n},从低到高sorting
  • mean_sin(v,x)= sum(sin(2 * pi * v_i / x),对于in {1,…,n})/ n
  • 在{1,…,n})/ n中,mean_cos(v,x)= sum(cos(2 * pi * v_i / x)
  • gcd_appeal (v,x)= 1-sqrt(mean_sin(v,x)^ 2 +(mean_cos(v,x)-1)^ 2)/ 2,其产生区间[0,1]中的数字。

目标是find最大化吸引力的x 。 这里是你的例子[2.468,3.700,6.1699]的( gcd_appeal )图,你可以发现最佳的GCD在x = 1.2337899957639993 在这里输入图像说明