algorithm - 2つのベクトル間の最適/スケール/シフトを見つける

Question

関数f（x）を表す2つのベクトルと、別のベクトルf（a x + b）、つまりf（x）のスケーリングおよびシフトされたバージョンがあります。最適なスケールとシフトファクターを見つけたいと思います。

*最良-最小二乗誤差、最尤法などを使用します。

何か案は？

例えば：

f1 = [0;0.450541598502498;0.0838213779969326;0.228976968716819;0.91333736150167;0.152378018969223;0.825816977489547;0.538342435260057;0.996134716626885;0.0781755287531837;0.442678269775446;0];
f2 = [-0.029171964726699;-0.0278570165494982;0.0331454732535324;0.187656956432487;0.358856370923984;0.449974662483267;0.391341738643094;0.244800719791534;0.111797007617227;0.0721767235173722;0.0854437239807415;0.143888234591602;0.251750993723227;0.478953530572365;0.748209818420035;0.908044924557262;0.811960826711455;0.512568916956487;0.22669198638799;0.168136111568694;0.365578085161896;0.644996661336714;0.823562159983554;0.792812945867018;0.656803251999341;0.545799498053254;0.587013303815021;0.777464637372241;0.962722388208354;0.980537136457874;0.734416947254272;0.375435649393553;0.106489547770962;0.0892376361668696;0.242467741982851;0.40610516900965;0.427497319032133;0.301874099075184;0.128396341665384;0.00246347624097456;-0.0322120242872125]

* f（x）は元に戻せない場合があることに注意してください...

ありがとう、

オハド

Answer 1

各 についてf(x)、 の絶対値を取得し、f(x)そのサポートに対する確率質量関数と見なすことができるように正規化します。E[x]の期待値と分散を計算しますVar[x]。次に、私たちはそれを持っています

  E[a x + b] = a E[x] + b
Var[a x + b] = a^2 Var[x]

上記の式と と の既知の値を使用して、E[x]とVar[x]を計算aしbます。例の値f1とf2例から値を取得すると、次の Octave スクリプトがこの手順を実行します。

% Octave script
% f1, f2 are defined as given in your example
f1 = [zeros(length(f2) - length(f1), 1); f1];

save_f1 = f1; save_f2 = f2;

f1 = abs( f1 ); f2 = abs( f2 );
f1 = f1 ./ sum( f1 ); f2 = f2 ./ sum( f2 );

mean = @(x)sum(((1:length(x))' .* x));
var = @(x)sum((((1:length(x))'-mean(x)).^2) .* x);

m1 = mean(f1); m2 = mean(f2);
v1 = var(f1); v2 = var(f2)

a = sqrt( v2 / v1 ); b = m2 - a * m1;

plot( a .* (1:length( save_f1 )) + b, save_f1, ...
      1:length( save_f2 ), save_f2 );
axis([0 length( save_f1 )];

そして、出力は

Answer 2

これは、単純で効果的ですが、おそらくやや単純なアプローチです。

最初に、両方の関数を使用して一般的なインターポレーターを作成してください。そうすれば、指定されたデータ ポイント間で両方の関数を評価できます。3 次スプライン補間器を使用しました。これは、提供された滑らかな関数のタイプに対して十分に一般的であるように思われるためです (追加のツールボックスは必要ありません)。

次に、多数のポイントでソース関数 ("元の") を評価します。この数値は、input として受け取るインライン関数のパラメーターとしても使用しますX。ここで、

X = [a b] 

(のようにax+b)。任意の入力に対してX、このインライン関数は計算します

1. 同じ x 位置にあるターゲット関数の関数値ですが、それぞれ および によってスケーリングおよびオフセットされaますb。

2. 結果の関数値と、以前に計算したソース関数の値との二乗差の合計。

このインライン関数はfminsearch、最初の見積もり (視覚的または自動手段で取得したもの) で使用します。あなたが提供した例では、いくつかのランダムなものを使用しましたが、それらはすべてほぼ最適な適合に収束しました.

上記のコードのすべて：

function s = findScaleOffset

    %% initialize

    f2 = [0;0.450541598502498;0.0838213779969326;0.228976968716819;0.91333736150167;0.152378018969223;0.825816977489547;0.538342435260057;0.996134716626885;0.0781755287531837;0.442678269775446;0];
    f1 = [-0.029171964726699;-0.0278570165494982;0.0331454732535324;0.187656956432487;0.358856370923984;0.449974662483267;0.391341738643094;0.244800719791534;0.111797007617227;0.0721767235173722;0.0854437239807415;0.143888234591602;0.251750993723227;0.478953530572365;0.748209818420035;0.908044924557262;0.811960826711455;0.512568916956487;0.22669198638799;0.168136111568694;0.365578085161896;0.644996661336714;0.823562159983554;0.792812945867018;0.656803251999341;0.545799498053254;0.587013303815021;0.777464637372241;0.962722388208354;0.980537136457874;0.734416947254272;0.375435649393553;0.106489547770962;0.0892376361668696;0.242467741982851;0.40610516900965;0.427497319032133;0.301874099075184;0.128396341665384;0.00246347624097456;-0.0322120242872125];

    figure(1), clf, hold on

    h(1) = subplot(2,1,1); hold on
    plot(f1);
    legend('Original')

    h(2) = subplot(2,1,2); hold on
    plot(f2);

    linkaxes(h)
    axis([0 max(length(f1),length(f2)), min(min(f1),min(f2)),max(max(f1),max(f2))])


    %% make cubic interpolators and test points

    pp1 = spline(1:numel(f1), f1);
    pp2 = spline(1:numel(f2), f2);

    maxX = max(numel(f1), numel(f2));
    N  = 100 * maxX;

    x2 = linspace(1, maxX, N);
    y1 = ppval(pp1, x2);

    %% search for parameters

    s = fminsearch(@(X) sum( (y1 - ppval(pp2,X(1)*x2+X(2))).^2 ), [0 0])

    %% plot results

    y2 = ppval( pp2, s(1)*x2+s(2));

    figure(1), hold on
    subplot(2,1,2), hold on    
    plot(x2,y2, 'r')
    legend('before', 'after')



end

結果：

s =
2.886234493867320e-001    3.734482822175923e-001

これは、データを生成したものとは逆の変換を計算することに注意してください。数字を逆にする：

>> 1/s(1) 
ans =    
    3.464721948700991e+000   % seems pretty decent 
>> -s(2)
ans = 
    -3.734482822175923e-001  % hmmm...rather different from 7/11!

(あなたが提供した 7/11 の値について確信が持てません。プロットを作成するために指定した正確な値を使用すると、ソース関数の近似値の精度が低下します...7/11 について確信がありますか?)

精度は次のいずれかで改善できます

1. 別のオプティマイザーを使用する ( fmincon、fminuncなど)
2. より高い精度が求めfminsearchられるoptimset
3. f1両方でより多くのサンプルポイントを持ちf2、補間の品質を向上させます
4. より良い初期見積もりの​​使用

とにかく、このアプローチはかなり一般的で、良い結果が得られます。また、ツールボックスも必要ありません。

ただし、これには大きな欠点が 1 つあります。見つかったソリューションがグローバル オプティマイザーではない可能性があります。たとえば、この方法の結果の品質は、提供する最初の見積もりに大きく左右される可能性があります。したがって、常に (差) プロットを作成して、最終的な解が正確であることを確認するか、そのような作業が多数ある場合は、何らかの品質係数を計算し、その上で最適化を別の方法で再開することにします。初期見積もり。

もちろん、フーリエ + メリン変換の結果 (以下の chaohuang が示唆するように) をこの方法の初期推定として使用することは非常に可能です。あなたが提供する単純な例ではやり過ぎかもしれませんが、これが実際に非常に役立つ状況を簡単に想像できます。

Answer 3

aサンプルデータで機能するスケールを推定するための非常に簡単なアプローチを次に示します。

a = length(f2) / length(f1)

これにより、3.4 という指定値に近い 3.4167 が得られます。その推定値が十分であれば、相関関係を使用してシフトを推定できます。

これはまさにあなたが求めたものではないことは承知していますが、データによっては許容できる代替手段である可能性があります。

Answer 4

Rody Oldenhuis と jstarr の両方の答えは正しいです。物事を要約し、それらを結び付けるためだけに、私自身の答えを追加しています。私はロディのコードを少し台無しにして、次のようになりました:

function findScaleShift
load f1f2

x0 = [length(f1)/length(f2) 0]; %initial guess, can do better
n=length(f1);
costFunc = @(z) sum((eval_f1(z,f2,n)-f1).^2);
opt.TolFun = eps; 
xopt=fminsearch(costFunc,x0,opt);
f1r=eval_f1(xopt,f2,n);
subplot(211);
plot(1:n,f1,1:n,f1r,'--','linewidth',5)
title(xopt);
subplot(212);
plot(1:n,(f1-f1r).^2);
title('squared error')
end


function y = eval_f1(x,f2,n)
t = maketform('affine',[x(1) 0 x(2); 0 1 0 ; 0 0 1]');
y=imtransform(f2',t,'cubic','xdata',[1 n ],'ydata',[1 1])';
end

これはゼロの結果をもたらします: この方法は正確ですが、網羅的であり、時間がかかる場合があります。もう 1 つの欠点は、極小値しか見つからず、最初の推定値 (x0) が遠い場合に誤った結果が返される可能性があることです。

一方、jstarr メソッドでは、次の結果が得られました。

xopt = [ 3.49655562549115         -0.676062367063033]

これは正解からの偏差が 10% です。かなり速い解決策ですが、私が要求したほど正確ではありませんが、それでも注意する必要があります。最良の結果を得るには、jstarr メソッドを Rody の目的とするメソッドの初期推測として使用して、正確な解を与える必要があると思います。

オハッド