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.

したがって、最初の 2 つの数値の疑似 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 になります。これらの値の小数部は、最小値が基本周波数にどれだけ近いかによって決定される N の値に対して、1/N に非常に近くなります。十分に洗練された一致が見つかるまで、N を増やしてループすることをお勧めします。この場合、N=1 (つまり、X=2.468 が基本周波数であると仮定) の場合、標準偏差は 0.3333 (3 つの値のうち 2 つが X * 1 から 0.5 ずれている) であり、許容できないほど高い値です。N=2 の場合 (つまり、2.468/2 が基本周波数であると仮定すると)、実質的にゼロの標準偏差が見つかります (3 つの値はすべて、X/2 の倍数の .001 以内です)。したがって、2.468/2 が近似値です。 GCD。

私の計画の主な欠点は、最小の測定値が最も正確な場合に最適に機能するということですが、そうではない可能性があります。これは、操作全体を複数回実行し、毎回測定値リストの最小値を破棄し、各パスの結果リストを使用してより正確な結果を決定することで軽減できます。結果を改善するもう 1 つの方法は、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 よりもはるかに大きな整数を使用していることに気付くかもしれません。 .

どの数字で割るかは大した問題ではありません。(たとえば、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 の Mathworldには、この問題の詳細な数学的処理がありますが、かなり簡単な説明は、グーグルで検索すると見つかります。

于 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 で開始し、間隔 I2 と I3 が整数を含むまで n1 をインクリメントします。これは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)

別の方法として、f が特定の周波数を下回る (n1 が上限に達する) まで、反復を継続し、min(J(n1))*n1 などの別の基準を最小化しようとすることもできます...

于 2012-11-15T00:27:36.513 に答える