Simulating estimators for the standard deviation

1

The standard deviation of x(:) is commonly computed as sqrt(sum((x-mean(x))**2)/(n-1)) where n = size(x), but Wikipedia says that instead of (n - 1), (n - 1.5) or (n - 1.5 + 1/(8*(n-1.0))) are more accurate. I have simulated this with random normal deviates for n = 100 but found that the usual (n - 1) formula (method 2 below) works best. For the function

pure function sd(x,imethod) result(xsd)
! formulas to compute standard deviation from Wikipedia https://en.wikipedia.org/wiki/Standard_deviation
real(kind=dp), intent(in) :: x(:)
integer      , intent(in) :: imethod
real(kind=dp)             :: xsd
real(kind=dp)             :: xmean,ssd
integer                   :: n
n = size(x)
if (n < 2) then
   xsd = -1.0_dp
   return
end if
xmean = mean(x)
ssd   = sum((x-xmean)**2)
select case (imethod)
   case (1)    ; xsd = sqrt(ssd/n)
   case (2)    ; xsd = sqrt(ssd/(n - 1))
   case (3)    ; xsd = sqrt(ssd/(n - 1.5_dp))
   case (4)    ; xsd = sqrt(ssd/(n - 1.5_dp + 1.0_dp/(8*n-8)))
   case default; xsd = -2.0_dp
end select
end function sd

with

main program 
program xstdev
use kind_mod    , only: dp
use ziggurat_mod, only: zigset, rnor
use stdev_mod   , only: mean,sd
implicit none
integer      , parameter :: ntotal = 10**7, ngroups = ntotal/100, n = ntotal/ngroups, &
                            nmethods = 4, niter = 10
real(kind=dp), parameter :: true_sd = 1.0_dp
logical      , parameter :: print_each = .false.
real(kind=dp)            :: x(ntotal),xsd(ngroups,nmethods),avg_err(niter,nmethods)
integer                  :: i,i1,i2,imethod,iter
print*,"#total, #groups, group_size =",ntotal,ngroups,n
call zigset(123456)
do iter=1,niter
   x = rnor(ntotal)
   x = x - mean(x)        ! set mean to zero
   x = x/sqrt(mean(x**2)) ! set variance to one
   i1 = 1
   do i=1,ngroups
      i2 = i1 + n - 1
      forall (imethod=1:nmethods) xsd(i,imethod) = sd(x(i1:i2),imethod) ! compute sd for subset using various methods
      i1 = i1 + n
   end do
   if (print_each) write (*,"(/,2a10)") "method","avg_error"
   do imethod=1,nmethods
      avg_err(iter,imethod) = mean(abs(xsd(:,imethod)-true_sd)) ! average absolute error of sd estimates for subsets
      if (print_each) write (*,"(i10,f10.6)") imethod,avg_err(iter,imethod)
   end do
end do
write (*,"(/,'OVERALL',/,2a10)") "method","avg_error"
do imethod=1,nmethods
   write (*,"(i10,f10.6)") imethod,mean(avg_err(:,imethod))
end do
end program xstdev

the average absolute error in the estimate of standard deviation is

OVERALL
    method avg_error
         1  0.056745
         2  0.056695
         3  0.056778
         4  0.056778

The full program is at GitHub. Have I coded the estimators of standard deviation incorrectly, or is the design of the simulation incorrect?

Edit: With

revised main program 
program xstdev
use kind_mod    , only: dp
use ziggurat_mod, only: zigset, rnor
use stdev_mod   , only: mean, rms, sd
implicit none
integer      , parameter :: ntotal = 10**7, ngroups = ntotal/100, n = ntotal/ngroups, &
                            nmethods = 4, niter = 10
real(kind=dp), parameter :: true_sd = 1.0_dp
logical      , parameter :: print_each = .false., standardize = .true.
real(kind=dp)            :: x(ntotal),xsd(ngroups,nmethods),avg_sd(niter,nmethods),avg_err(niter,nmethods),rmse(niter,nmethods)
integer                  :: i,i1,i2,imethod,iter
print*,"#total, #groups, group_size =",ntotal,ngroups,n
call zigset(123456)
do iter=1,niter
   x = rnor(ntotal)
   if (standardize) then
      x = x - mean(x)        ! set mean to zero
      x = x/sqrt(mean(x**2)) ! set variance to one
   end if
   i1 = 1
   do i=1,ngroups
      i2 = i1 + n - 1
      forall (imethod=1:nmethods) xsd(i,imethod) = sd(x(i1:i2),imethod) ! compute sd for subset using various methods
      i1 = i1 + n
   end do
   if (print_each) write (*,"(/,*(a14))") "method","avg_sd","rmse","avg_error"
   do imethod=1,nmethods
      avg_sd(iter,imethod)  = mean(xsd(:,imethod))
      avg_err(iter,imethod) = mean(abs(xsd(:,imethod)-true_sd)) ! average absolute error of sd estimates for subsets
      rmse(iter,imethod)    = rms(xsd(:,imethod)-true_sd)
      if (print_each) write (*,"(i14,*(f14.10))") imethod,avg_sd(iter,imethod),rmse(iter,imethod),avg_err(iter,imethod)
   end do
end do
write (*,"(/,'OVERALL',/,4a14)") "method","avg_sd","rmse","avg_error"
do imethod=1,nmethods
   write (*,"(i14,*(f14.10))") imethod,mean(avg_sd(:,imethod)),mean(rmse(:,imethod)),mean(avg_err(:,imethod))
end do
end program xstdev

I do find that method 3 gives the least biased estimate of the standard deviation, although estimators in 1 and 2 still have a lower RMSE:

 #total, #groups, group_size =    10000000      100000         100

OVERALL
        method        avg_sd          rmse     avg_error
             1  0.9924683782  0.0002245609  0.0567451915
             2  0.9974682505  0.0002245620  0.0566951266
             3  0.9999966912  0.0002249881  0.0567783850
             4  0.9999902820  0.0002249867  0.0567780829

Since method 2 with the (n-1) denominator is common, I think I will continue to use it.

2

I have never seen alternatives like that - it reminds me of the legion “plotting positions” for QQ-plots and the like. But IIRC the (n-1) factor is due to Bessel, not the least name.

3

See the discussion around equation 14.1.8 in Press et al, rounding errors and how to mitigate

4

Thanks for the suggestion. I tried adding the two-pass correction term to get

ssd = sum((x-xmean)**2) - ((sum(x-xmean))**2)/n

but the results are unchanged, so I don’t think the results are due to rounding errors.

5

The results seem consistent with the Wikipedia page (under Estimation). Compared with the conventional Method 2 (division by (n-1)), we get that.

  • Method 1 (division by n) has lower RMSE
  • Methods 3 and 4 have less bias.

I will also keep using Method 2 in general, but it’s a good reminder that if you want to optimize some particular property of the errors, there can be better choices.

1 Like
6

Thanks, I should have read the full article. It appears that division by (n - 0.5) (denoted Method 5 below) has both lower RMSE and bias than division by n. The simulation results are


OVERALL
        method        avg_sd          rmse     avg_error
             1  0.9924683782  0.0002245609  0.0567451915
             2  0.9974682505  0.0002245620  0.0566951266
             3  0.9999966912  0.0002249881  0.0567783850
             4  0.9999902820  0.0002249867  0.0567780829
             5  0.9949588924  0.0002244215  0.0566846522

I tried (n-c) for c = 0.1, 0.2, ..., 0.9, and c = 0.5 still gave the lowest RMSE. It also does so for groups of 10 rather than 100 observations. The revised program is at GitHub. A problem with using (n-0.5), besides its lack of theoretical foundation (other than being the midpoint of two common methods), is that one’s results results would not be exactly replicable using the the built-in standard deviation functions of R, NumPy, Matlab etc.