From c8218c1bc42732cd09e422d27123719ffe027b2d Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sun, 11 Jul 2021 06:22:44 +0200 Subject: [PATCH] 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. --- .../@MOTSimulator/bootstrapErrorEstimation.m | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m diff --git a/MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m b/MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m new file mode 100644 index 0000000..731f8b4 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m @@ -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 \ No newline at end of file