From d050c72d1d7e12df6fdc4d2ede24f481950f2265 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Mar 2012 14:16:29 -0400 Subject: [PATCH] Added Bessel's correction to make estimate of variance unbiased. --- src/intercycle.F90 | 7 +++++-- src/tally.F90 | 8 ++++++-- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 97a7b9f86d..3d5f0cd704 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -394,12 +394,15 @@ contains ! Accumulate single batch realizations of k call accumulate_batch_estimate(global_tallies) - ! Determine mean and standard deviation of mean + ! Determine sample mean of keff keff = global_tallies(K_ANALOG) % sum/n - keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n) ! Display output for this batch if (current_batch > n_inactive + 1) then + ! Determine standard deviation of the sample mean of keff + keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - & + keff*keff)/(n - 1)) + if (entropy_on) then write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, & entropy, keff, keff_std diff --git a/src/tally.F90 b/src/tally.F90 index 555e873ad2..4702ed9baf 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1532,9 +1532,13 @@ contains type(TallyScore), intent(inout) :: score - ! Calculate mean and standard deviation of the mean + ! Calculate sample mean and standard deviation of the mean -- note that we + ! have used Bessel's correction so that the estimator of the variance of the + ! sample mean is unbiased. + score % sum = score % sum/n_active - score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2)/n_active) + score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & + (n_active - 1)) end subroutine calculate_statistics