Major change how the loading rate is calculated from the raw output of the simulator - Zero crossing of the autocorrelation of the time series - which is the number of loaded atoms at each time step - is used to determine the length of each of a 1000 bootstrap samples obtained by sampling with replacement. Mean in each sample is calculated and a sampling distribution of means obtained. The capture ratio is the mean of this distribution which when multiplied by the "reduced" flux gives the loading rate. Standard Error and 95% confidence interval are also calculated from the sampling distribution.

This commit is contained in:
Karthik 2021-07-11 06:22:44 +02:00
parent 1a8975c218
commit c8218c1bc4

View File

@ -0,0 +1,40 @@
function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, 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;
MeanLoadingRatioInEachSample = zeros(1,NumberOfBootsrapSamples);
for SampleNumber = 1:NumberOfBootsrapSamples
BoostrapSample = datasample(NumberOfLoadedAtoms, SampleLength); % Sample with replacement
MeanLoadingRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
end
LoadingRate = mean(MeanLoadingRatioInEachSample) * this.ReducedFlux;
Variance = 0; % Bootstrap Estimate of Variance
for SampleNumber = 1:NumberOfBootsrapSamples
Variance = Variance + (MeanLoadingRatioInEachSample(SampleNumber) - mean(MeanLoadingRatioInEachSample))^2;
end
StandardError = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * this.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