找到两个向量之间的最佳/比例/转换

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]和Var[x]来计算a和b.以你的价值观f1,并f2从你的榜样,下面的八度脚本执行此过程:

% 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所示)作为该方法的初始估计.对于您提供的简单示例,这可能有点过分,但我可以很容易地想象出这可能非常有用的情况.