From 792eb39e89d711640efd6b199e0cbe359d3bcb9c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 20 Aug 2012 18:18:00 -0400 Subject: [PATCH] Fixed degrees of freedom for confidence intervals using t_percentile. --- src/intercycle.F90 | 4 ++-- src/output.F90 | 2 +- src/state_point.F90 | 2 +- src/tally.F90 | 2 +- src/utils/statepoint_histogram.py | 2 +- src/utils/statepoint_meshplot.py | 4 ++-- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index b2729df997..16a61f015c 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -407,7 +407,7 @@ contains if (confidence_intervals) then ! Calculate t-value for confidence intervals alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations) + t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) else t_value = ONE end if @@ -441,7 +441,7 @@ contains if (confidence_intervals) then ! Calculate t-value for confidence intervals alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations) + t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) else t_value = ONE end if diff --git a/src/output.F90 b/src/output.F90 index 97ec86cfe4..dc4df8734c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1237,7 +1237,7 @@ contains if (confidence_intervals) then ! Calculate t-value for confidence intervals alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations) + t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) ! Adjust sum_sq global_tallies(:) % sum_sq = t_value * global_tallies(:) % sum_sq diff --git a/src/state_point.F90 b/src/state_point.F90 index bd400faf63..90f7ee5418 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -347,7 +347,7 @@ contains if (confidence_intervals) then ! Calculate t-value for confidence intervals alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations) + t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) else t_value = ONE end if diff --git a/src/tally.F90 b/src/tally.F90 index 92f9c5b1f7..f7ac74beec 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1864,7 +1864,7 @@ contains ! Calculate t-value for confidence intervals if (confidence_intervals) then alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations) + t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) end if do i = 1, n_tallies diff --git a/src/utils/statepoint_histogram.py b/src/utils/statepoint_histogram.py index da86e9cecd..c1144d7251 100755 --- a/src/utils/statepoint_histogram.py +++ b/src/utils/statepoint_histogram.py @@ -21,7 +21,7 @@ sp = StatePoint(filename) # Calculate t-value for 95% two-sided CI n = sp.current_batch - sp.n_inactive -t_value = scipy.stats.t.ppf(0.975, n) +t_value = scipy.stats.t.ppf(0.975, n - 1) # Loop over all tallies for i, t in enumerate(sp.tallies): diff --git a/src/utils/statepoint_meshplot.py b/src/utils/statepoint_meshplot.py index 0039ff6d73..0b2278c96c 100755 --- a/src/utils/statepoint_meshplot.py +++ b/src/utils/statepoint_meshplot.py @@ -22,7 +22,7 @@ sp = StatePoint(filename) # Calculate t-value for 95% two-sided CI n = sp.current_batch - sp.n_inactive -t_value = scipy.stats.t.ppf(0.975, n) +t_value = scipy.stats.t.ppf(0.975, n - 1) # Loop over all tallies print("Reading data...") @@ -60,7 +60,7 @@ for t in sp.tallies: # Make figure print("Making colorplot...") - plt.imshow(error, interpolation="nearest") + plt.imshow(mean, interpolation="nearest") plt.colorbar() plt.xlim((0,nx)) plt.ylim((0,ny))