MATLAB 中 FOR 循环中的回归

Question

我有以下代码：

colBIN = {0.050, 0.055, 0.060, 0.065, 0.070, 0.075, 0.080, 0.085, 0.090, 0.095,0.1};

for i = 1 : length(colBIN)-1
    colBIN{i,2} = find(cols(:,1) <= cell2mat(colBIN(i+1,1)) & cols(:,1) > cell2mat(colBIN(i,1)));
end

rowBIN = {0.045, 0.046, 0.047, 0.048, 0.049, 0.050, 0.051, 0.052};

for i = 1 : length(rowBIN)-1
    rowBIN{i,2} = find(rows(:,1) <= cell2mat(rowBIN(i+1,1)) & rows(:,1) > cell2mat(rowBIN(i,1))); 
end

binCombos = cell(length(rowBIN)-1,length(colBIN)-1);

for m = 1 : length(rowBIN)-1
    for n = 1 : length(colBIN)-1
        binCombos{n,m} = intersect( rowBIN{m,2}(:,1),colBIN{n,2}(:,1));
    end
end


binRows = size(binCombos,1);
binCols = size(binCombos,2)-1;

j = j + 1;
for n = 1 : binRows; 
    for m = 1 : binCols;
       thisBin = binCombos{n,m}(:,:); 
       if isempty(thisBin)==0

       %polyfit
       quadmod = polyfit(x_vrbl(thisBin), y_vrbl(thisBin), 2);
       interval = 0.0:0.001:1;
       quadmodcurve = polyval(quadmod,interval); 
       [r2 rmse] = rsquare(y_vrbl(thisBin), quadmodcurve); 
       plot(x_vrbl(thisBin), y_vrbl(thisBin), '*', interval, quadmodcurve);
       xlabel('x_vrbl');
       ylabel('y_vrbl');
       axis([0,1,0,1]);
       header = ['R^2 =' num2str(r2),'coeffs:',num2str(quadmod)];
       title(header);
       saveas(gcf, sprintf('plot_%d.pdf', j));

       %residuals
       res = y_vrbl(thisBin) - quadmodcurve;
       plot(x_vrbl(thisBin),res,'+');
       header2 = ['residuals'];
       title(header2);
       saveas(gcf, sprintf('residuals_%d.pdf', j));

       end
       j = j + 1;
   end
end

说明/问题：

binCombos 是一个二维元胞数组，每个元胞具有不均匀数量的数据点。我正在为每个唯一单元格的数据拟合一条二次曲线，并尝试（不成功）输出 R^2 值以及 绘制残差。

我认为问题与以下事实有关：polyval 函数所需的“间隔”在尝试查找 rsquare 时与 y_vrbl(thisBin) 的数组大小不匹配，同样用于计算残差。例如，如果我设置了interval = x_vrbl(thisBin)，那么残差“起作用”，但 polyfit 全部搞砸了。

Answer 1

我的猜测是这应该可行：

quadmodcurve = polyval(quadmod,y_vrbl(thisBin)); 
[r2 rmse] = rsquare(y_vrbl(thisBin), quadmodcurve);
interval = 0.0:0.001:1;
quadmodcurve = polyval(quadmod,interval); 

为了确定拟合质量，您必须仅在样本的 x 值处评估多项式。为了绘制完整的多项式图，您需要在更多且规则间隔的 x 值上对其进行评估。

Answer 2

我设法让您的代码使用http://dropproxy.com/f/4B6 的数据和file exchange 的 rsquare 函数运行 纠正一些错误后：

d = importdata('sample_data.xlsx');
y_vrbl = d.data(:, 1);
x_vrbl = d.data(:, 2);
rows = d.data(:, 3);
cols = d.data(:, 4);

cb = {0.050, 0.055, 0.060, 0.065, 0.070, 0.075, 0.080, 0.085, 0.090, 0.095,0.1};

for i = 1 : length(cb)-1
    colBIN{i,2} = find(cols(:,1) <= cell2mat(cb(i+1)) & cols(:,1) > cell2mat(cb(i)));
end

rb = {0.045, 0.046, 0.047, 0.048, 0.049, 0.050, 0.051, 0.052};

for i = 1 : length(rb)-1
    rowBIN{i,2} = find(rows(:,1) <= cell2mat(rb(i+1)) & rows(:,1) > cell2mat(rb(i)));
end

binCombos = cell(length(rowBIN)-1,length(colBIN)-1);

for m = 1 : length(rowBIN)-1
    for n = 1 : length(colBIN)-1
        binCombos{n,m} = intersect( rowBIN{m,2}(:,1),colBIN{n,2}(:,1));
    end
end


binRows = size(binCombos,1);
binCols = size(binCombos,2)-1;

j = 1;
for n = 1 : binRows;
    for m = 1 : binCols;
        thisBin = binCombos{n,m}(:,:);
        if ~isempty(thisBin)

            % polyfit
            quadmod = polyfit(x_vrbl(thisBin), y_vrbl(thisBin), 2);

            % compute residuals and R²
            quadmodcurve = polyval(quadmod,y_vrbl(thisBin));
            [r2, rmse] = rsquare(y_vrbl(thisBin), quadmodcurve);
            res = y_vrbl(thisBin) - quadmodcurve;

            % plot fit
            interval = 0.0:0.001:1;
            quadmodcurve = polyval(quadmod,interval);
            plot(x_vrbl(thisBin), y_vrbl(thisBin), '*', interval, quadmodcurve);
            xlabel('x_vrbl');
            ylabel('y_vrbl');
            axis([0,1,0,1]);
            header = ['R^2 =' num2str(r2),'coeffs:',num2str(quadmod)];
            title(header);
            saveas(gcf, sprintf('plot_%d.pdf', j));

            % plot residuals
            plot(x_vrbl(thisBin),res,'+');
            header2 = ['residuals'];
            title(header2);
            saveas(gcf, sprintf('residuals_%d.pdf', j));

        end
        j = j + 1;
    end
end

对我来说，拟合看起来不错，但在大多数情况下，线性函数可能就足够了，而二次项不是必需的。

关于您剩下的问题：我不是使用 R² 进行非线性拟合的专家（请参阅 Coefficient of determination 上的注释 2），但您使用的实现对我来说似乎有点可疑。大多数情况下输出 0 的原因是 rsquare.m 的第 65 行上的 max 函数，它可以防止返回负值。由于您的多项式拟合确实包含一个常数项，因此将该函数称为

[r2, rmse] = rsquare(y_vrbl(thisBin), quadmodcurve, false);

似乎更合适，在大多数情况下 R² > 0.9。

我的建议：检查 R² 是否是您案例中拟合优度的正确度量，并检查该函数是否正确实现它。 Matlab 附带的函数可以开箱即用地被信任，但 Matlab 文件交换上的帖子没有质量保证。