algorithm - 找到两个向量之间的最佳/缩放/偏移

Question

我有两个向量表示函数 f(x)，另一个向量 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]和。从您的示例中获取您的值，以下 Octave 脚本执行此过程：Var[x]abf1f2

% 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

这是一种简单、有效但可能有点幼稚的方法。

首先确保通过这两个函数创建一个通用插值器。这样您就可以在给定的数据点之间评估这两个函数。我使用了三次样条插值器，因为这对于您提供的平滑函数类型来说似乎足够通用（并且不需要额外的工具箱）。

然后在大量点上评估源函数（“原始”）。将此数字也用作内联函数中的参数，该函数作为输入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. 使用更好的初始估计

无论如何，这种方法非常通用，并给出了很好的结果。它也不需要工具箱。

但它有一个主要缺点——找到的解决方案可能不是全局优化器，例如，这种方法的结果质量可能对您提供的初始估计非常敏感。因此，请始终制作（差异）图以确保最终解决方案是准确的，或者如果您有大量此类事情要做，请计算某种品质因数，您决定在此基础上使用不同的方法重新开始优化初步估计。

当然，很可能使用傅里叶+梅林变换的结果（如下面的 chaohuang 所建议的）作为对该方法的初始估计。对于您提供的简单示例来说，这可能有点过头了，但我可以很容易地想象出这确实非常有用的情况。

Answer 3

对于比例因子 a，您可以通过计算两个信号的幅度谱的比率来估计它，因为傅里叶变换是不变的。

类似地，您可以通过使用尺度不变的梅林变换来估计移位因子 b。

Answer 4

这是一种超级简单的方法来估计适用于您的示例数据的规模a：

a = length(f2) / length(f1)

这给出了 3.4167，接近于您声明的 3.4 值。如果该估计值足够好，您可以使用相关性来估计偏移。

我意识到这并不完全是您所要求的，但根据数据，它可能是一个可以接受的替代方案。

Answer 5

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

结果为零： 此方法准确但详尽，可能需要一些时间。另一个缺点是它只能找到局部最小值，如果初始猜测 (x0) 很远，可能会给出错误的结果。

另一方面，jstarr 方法给出了以下结果：

xopt = [ 3.49655562549115         -0.676062367063033]

与正确答案有 10% 的偏差。相当快的解决方案，但不如我要求的准确，但仍应注意。我认为为了获得最佳结果，应该使用 jstarr 方法作为 Rody 所用方法的初始猜测，给出准确的解决方案。

奥哈德