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