function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, ovenObj,  NumberOfLoadedAtoms)
    n                  = this.NumberOfAtoms;
    NumberOfTimeSteps  = int64(this.SimulationTime/this.TimeStep);
    
    Autocorrelation = autocorr(NumberOfLoadedAtoms,'NumLags', double(NumberOfTimeSteps - 1));

    if Autocorrelation(1)~=0 
        CorrelationFactor                  = table(Helper.findAllZeroCrossings(linspace(1, double(NumberOfTimeSteps), double(NumberOfTimeSteps)), Autocorrelation)).Var1(1);
        if ~isnan(CorrelationFactor)
            SampleLength                   = floor(CorrelationFactor);
            NumberOfBootsrapSamples            = 1000;
            MeanCaptureRatioInEachSample       = zeros(1,NumberOfBootsrapSamples);
            for SampleNumber = 1:NumberOfBootsrapSamples
                BoostrapSample                             = datasample(NumberOfLoadedAtoms, SampleLength); % Sample with replacement
                MeanCaptureRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
            end

            LoadingRate = mean(MeanCaptureRatioInEachSample) * ovenObj.ReducedFlux;

            Variance = 0; % Bootstrap Estimate of Variance
            for SampleNumber = 1:NumberOfBootsrapSamples
                Variance = Variance + (MeanCaptureRatioInEachSample(SampleNumber) - mean(MeanCaptureRatioInEachSample))^2;
            end

            StandardError      = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * ovenObj.ReducedFlux;

            ts                 = tinv([0.025  0.975],NumberOfBootsrapSamples-1);  % T-Score
            ConfidenceInterval = LoadingRate + ts*StandardError;                  % 95% Confidence Intervals
        else
            LoadingRate   = nan;
            StandardError = nan;
            ConfidenceInterval = [nan nan];
        end
    else
        LoadingRate   = nan;
        StandardError = nan;
        ConfidenceInterval = [nan nan];
    end
end