How to optimize my code? I'm generating resampling distributions of estimates obtained from fitglm function?

I want to optimize my code. I want to generate a distribution (n=1000) of Estimates form a GLM in which I shuffle or use a cicle shift in my "y" variable, to compare it with my real data.

I'm using fitglm to predict y. In my case I have a matrix of 5 variables predictors (2736x5 double) and y (2736x1 double).

I want to know if there is possible to use bootsrp function in my code or how can I optimize it, instead of using a for loop.

% My data 
mdl = fitglm(predictors,y,'linear', ...
    'Distribution','poisson','PredictorVars',varNames,...
    'CategoricalVars',[1,2,3,4]);


nIter = 1000;
minShift = 2;
maxShift = round(length(y)/2);


for k = 1 : nIter
    %Shuffling
    yShuffled = y(randperm(length(y)));
    shuffledmdl = fitglm(predictors,yShuffled,'linear', ...
        'Distribution','poisson','PredictorVars',varNames,...
        'CategoricalVars',[1,2,3,4]);
    
    coeffValShuffle(k,1) = shuffledmdl.Coefficients.Estimate(1);
    coeffValShuffle(k,2) = shuffledmdl.Coefficients.Estimate(2);
    coeffValShuffle(k,3) = shuffledmdl.Coefficients.Estimate(3);
    coeffValShuffle(k,4) = shuffledmdl.Coefficients.Estimate(4);
    coeffValShuffle(k,5) = shuffledmdl.Coefficients.Estimate(5);
    coeffValShuffle(k,6) = shuffledmdl.Coefficients.Estimate(6);
    coeffValShuffle(k,7) = shuffledmdl.Coefficients.Estimate(7);

    %Circle shift
    shifting = randi([minShift maxShift]);
    yShifted = circshift(y, shifting);
    circShiftmdl = fitglm(predictors,yShifted,'linear', ...
        'Distribution','poisson','PredictorVars',varNames,...
        'CategoricalVars',[1,2,3,4]);
    
    coeffValCirleShift(k,1) = circShiftmdl.Coefficients.Estimate(1);
    coeffValCirleShift(k,2) = circShiftmdl.Coefficients.Estimate(2);
    coeffValCirleShift(k,3) = circShiftmdl.Coefficients.Estimate(3);
    coeffValCirleShift(k,4) = circShiftmdl.Coefficients.Estimate(4);
    coeffValCirleShift(k,5) = circShiftmdl.Coefficients.Estimate(5);
    coeffValCirleShift(k,6) = circShiftmdl.Coefficients.Estimate(6);
    coeffValCirleShift(k,7) = circShiftmdl.Coefficients.Estimate(7);
end

to predict y.

Answer 1

[Karan Singh]

Hi @María ML,

I dont seee a way to directly integrate the "bootstrp" function into the code. However you can use a few imporvements.

1. The loop can be parallelized using "parfor" to speed things up. https://in.mathworks.com/help/parallel-computing/parfor.html
2. MATLAB is inefficient with dynamically growing arrays inside loops. Preallocating would avoids costly memory reallocation.
3. Instead of manually specifying "coeffValShuffle(k,1:7)", use "numCoeffs = size(mdl.Coefficients.Estimate, 1)". So that you dont reach ut for the same again and again.

% My data 
mdl = fitglm(predictors, y, 'linear', ...
    'Distribution', 'poisson', 'PredictorVars', varNames, ...
    'CategoricalVars', [1,2,3,4]);

nIter = 1000;
minShift = 2;
maxShift = round(length(y) / 2);
numCoeffs = size(mdl.Coefficients.Estimate, 1);  % Dynamically get coefficient count

% Preallocate matrices for efficiency
coeffValShuffle = zeros(nIter, numCoeffs);
coeffValCircleShift = zeros(nIter, numCoeffs);

% Use parallel computing for speedup
parpool; % Start parallel pool if not already started
parfor k = 1:nIter
    % Shuffling
    yShuffled = y(randperm(length(y)));
    shuffledmdl = fitglm(predictors, yShuffled, 'linear', ...
        'Distribution', 'poisson', 'PredictorVars', varNames, ...
        'CategoricalVars', [1,2,3,4]);
    
    coeffValShuffle(k, :) = shuffledmdl.Coefficients.Estimate';

    % Circular shift
    shifting = randi([minShift, maxShift]);
    yShifted = circshift(y, shifting);
    circShiftmdl = fitglm(predictors, yShifted, 'linear', ...
        'Distribution', 'poisson', 'PredictorVars', varNames, ...
        'CategoricalVars', [1,2,3,4]);
    
    coeffValCircleShift(k, :) = circShiftmdl.Coefficients.Estimate';
end
delete(gcp); % Close parallel pool

Karan

[María ML]

Thanks @Karan Singh, that would be very useful for me. Thanks for all the optimization comments.

Let me try this parfor funtion,