«Приближенный» наибольший общий делитель

Question

Предположим, у вас есть список чисел с плавающей запятой, которые приблизительно кратны общему количеству, например

2,468, 3,700, 6,1699

, которые примерно все кратны 1,234. Как бы вы охарактеризовали этот «приблизительный gcd», и как бы вы приступили к его вычислению или оценке?

Строго связан с моим ответом на этот вопрос .

Answer 1

Вы можете запустить алгоритм Евклида 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, а mutiples - 2,3,5. Чтобы улучшить этот результат, вы можете использовать линейную регрессию для точек данных:

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

Наклон улучшенного псевдо-gcd.

Предупреждение: первая часть этого алгоритма является численно нестабильной - если вы начинаете с очень грязных данных, у вас проблемы.

Answer 2

Выразите ваши измерения как кратные младшего. Таким образом, ваш список становится 1.00000, 1.49919, 2.49996. Дробные части этих значений будут очень близки к 1 / Nths, для некоторого значения N продиктовано тем, насколько близко ваше самое низкое значение к основной частоте. Я бы посоветовал циклически увеличивать N, пока вы не найдете достаточно точное соответствие. В этом случае для N = 1 (то есть, если X = 2,468 - ваша основная частота), вы найдете стандартное отклонение 0,3333 (два из трех значений - 0,5 от X * 1), что недопустимо высоко. Для N = 2 (то есть при условии, что 2,468 / 2 является вашей основной частотой) вы найдете стандартное отклонение практически равным нулю (все три значения находятся в пределах .001 от кратного X / 2), поэтому 2,468 / 2 - это ваше приблизительное значение. НОД.

Главный недостаток моего плана заключается в том, что он работает лучше всего, когда самое низкое измерение является наиболее точным, что, скорее всего, не так. Этого можно избежать, выполнив всю операцию несколько раз, каждый раз отбрасывая самое низкое значение в списке измерений, а затем используйте список результатов каждого прохода, чтобы определить более точный результат. Другим способом уточнения результатов будет настройка GCD, чтобы минимизировать стандартное отклонение между целочисленными коэффициентами GCD и измеренными значениями.

Answer 3

Это напоминает мне о проблеме нахождения хороших приближений рациональных чисел к действительным числам. Стандартная методика - расширение непрерывной дроби:

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

Мы могли бы применить это непосредственно к подходу Джонатана Леффлера и Спарра:

>>> 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. Все побеждают!

Answer 4

Я предполагаю, что все ваши числа кратны целому числу значениям. В остальной части моего объяснения 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, вы можете использовать линейную регрессию, чтобы минимизировать ошибку линейной модели теперь, когда у вас есть соответствующие зависимые переменные. Самый простой способ сделать это - использовать метод наименьших квадратов. Mathworld Вольфрама содержит углубленную математическую трактовку этой проблемы, но довольно простое объяснение можно найти с некоторым поиском в Google.

Answer 5

Интересный вопрос ... не просто.

Полагаю, я бы посмотрел на соотношение значений выборки:

• 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.

Answer 6

Решение, которое я видел и использовал сам, состоит в том, чтобы выбрать некоторую константу, скажем 1000, умножить все числа на эту константу, округлить их до целых чисел, найти GCD этих целых чисел, используя стандартный алгоритм, а затем разделить результат указанная константа (1000). Чем больше константа, тем выше точность.

Answer 7

Я нашел этот вопрос в поисках ответов на мой вопрос в MathStackExchange ( здесь и здесь ).

Мне удалось (пока) измерить апелляцию фундаментальной частоты с помощью списка частот гармоник (следуя номенклатуре звука / музыки), что может быть полезно, если у вас уменьшенное число вариантов и возможно рассчитать привлекательность каждого из них, а затем выбрать наиболее подходящий вариант.

C & P из моего вопроса в MSE (там форматирование красивее):

• будучи v списком {v_1, v_2, ..., v_n}, упорядоченным снизу вверх
• mean_sin (v, x) = сумма (sin (2 * pi * v_i / x), для i в {1, ..., n}) / n
• mean_cos (v, x) = сумма (cos (2 * pi * v_i / x), для i в {1, ..., n}) / n
• gcd_appeal (v, x) = 1 - sqrt (mean_sin (v, x) ^ 2 + (mean_cos (v, x) - 1) ^ 2) / 2, что дает число в интервал [0,1].

Цель состоит в том, чтобы найти х, который максимизирует апелляцию . Вот график ( gcd_appeal ) для вашего примера [2.468, 3.700, 6.1699], где вы обнаружите, что оптимальный GCD равен x = 1.2337899957639993

Answer 8

Это переформулировка решения shsmurfy, когда вы априори выбираете 3 положительных допуска (e1, e2, e3)
Задача состоит в том, чтобы искать наименьшие натуральные числа (n1, n2, n3) и, таким образом, наибольшую корневую частоту f, такую ​​что:

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

Мы предполагаем, что 0 <= f1 <= f2 <= f3 <br> Если мы исправим 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
Затем мы выбираем наименьшее целое число n2 в интервале I2 и наименьшее целое число n3 в интервале I3.

При нормальном распределении ошибок с плавающей запятой наиболее вероятной оценкой корневой частоты 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 не достигнет верхнего предела) ...