向量化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.