crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Fixed degrees of freedom for confidence intervals using t_percentile.

Browse source
This commit is contained in:
Paul Romano 2012-08-20 18:18:00 -04:00
parent a258b7f73f
commit 792eb39e89
6 changed files with 8 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/intercycle.F90
View file

@ -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

2
src/output.F90
View file

@ -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

2
src/state_point.F90
View file

@ -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

2
src/tally.F90
View file

@ -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

2
src/utils/statepoint_histogram.py
View file

@ -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):

4
src/utils/statepoint_meshplot.py
View file

@ -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))