向量化2d卡方網格搜索

Question

最佳擬合線性參數A和B（y = Ax + b）對應於這些參數上的卡方函數的最小值。我想做一個暴力網格搜索的全局最小值（保證，因爲2參數線性卡方是一個拋物面），並已實現它與3嵌套循環（下），但要避免循環（即矢量化使用數組屬性）。

卡方（加權最小二乘）被定義爲（僞代碼）：

卡方（K，J）=總和（值Y [i] - （A [k]的* X [I] + B [j]））/ yerr [I]）^ 2。

下面是Matlab代碼，填充100 x 100網格，參數值爲A和B（每個值爲100）的10,000個組合的卡方值。有三個數據陣列：x,y和yerr。

感謝您對兩參數線性卡方網格的空洞版本的任何幫助！

基思

% create parameter grid 
    a = linspace(85,110,100); 
    b = linspace(10,35,100); 
    [A,B] = meshgrid(a,b); 

    % calculate chi-square over parameter grid 
    chi2(100,100) = zeros; 

    for k = 1:100; 
     for j = 1:100; 
      for i = 1:length(y) 
      chi2a = ((y(i)-a(k)*x(i)-b(j))/yerr(i)).^2; 
      chi2(k,j) = chi2(k,j)+chi2a; 
      end 
     end 
    end 

Answer 1

我們可以bsxfun它 -

x3d = reshape(x,1,1,numel(x)); 
y3d = reshape(y,1,1,numel(y)); 
yerr3d = reshape(yerr,1,1,numel(yerr)); 
p0 = bsxfun(@minus, bsxfun(@minus,y3d,bsxfun(@times,a(:),x3d)), b); 
p1 = bsxfun(@rdivide, p0, yerr3d); 
out = sum(p1.^2,3); 

隨着MATLAB的隱式擴張，計算p0和p1將簡化爲 -

p0 = ((y3d - a(:).*x3d) - b); 
p1 = p0 ./yerr3d; 

計時 -

% Setup 
N = 2000; 
x = rand(N,1); 
y = rand(N,1); 
yerr = rand(N,1); 

a = linspace(85,110,100); 
b = linspace(10,35,100); 

我們得到 -

----------- With loopy method ------------------------- 
Elapsed time is 1.056787 seconds. 
----------- With BSXFUN method ------------------------- 
Elapsed time is 0.109601 seconds.