Jump to content

Cumulative standard deviation: Difference between revisions

(→‎Old style, four ways: Agh! More mistypes! And cut&pasted too.)
Line 1,118:
 
=={{header|Fortran}}==
{{works with|Fortran|2003 and later}}
{{incorrect|Fortran|It does not return the ''running'' standard deviation of the series.}}
<lang fortran>
{{trans|C}}
program standard_deviation
implicit none
integer(kind=4), parameter :: dp = kind(0.0d0)
 
real(kind=dp), dimension(:), allocatable :: vals
{{works with|Fortran|95 and later}}
integer(kind=4) :: i
 
real(kind=dp), dimension(8) :: sample_data = (/ 2, 4, 4, 4, 5, 5, 7, 9 /)
This one imitates C and suffers the same problems: the function is not thread-safe and must be used to compute the stddev for one set per time.
 
do i = lbound(sample_data, 1), ubound(sample_data, 1)
<lang fortran>program Test_Stddev
call sample_add(vals, sample_data(i))
implicit none
write(*, fmt='(''#'',I1,1X,''value = '',F3.1,1X,''stddev ='',1X,F10.8)') &
 
i, sample_data(i), stddev(vals)
real, dimension(8) :: v = (/ 2,4,4,4,5,5,7,9 /)
integer :: i
real :: sd
 
do i = 1, size(v)
sd = stat_object(v(i))
end do
 
if (allocated(vals)) deallocate(vals)
print *, "std dev = ", sd
 
contains
! Adds value :val: to array :population: dynamically resizing array
subroutine sample_add(population, val)
real(kind=dp), dimension(:), allocatable, intent (inout) :: population
real(kind=dp), intent (in) :: val
 
real(kind=dp), dimension(:), allocatable :: tmp
recursive function stat_object(a, cmd) result(stddev)
realinteger(kind=4) :: stddevn
real, intent(in) :: a
character(len=*), intent(in), optional :: cmd
 
if (.not. allocated(population)) then
real, save :: summa = 0.0, summa2 = 0.0
allocate(population(1))
integer, save :: num = 0
population(1) = val
 
real :: m
 
if ( .not. present(cmd) ) then
num = num + 1
summa = summa + a
summa2 = summa2 + a*a
stddev = stat_object(0.0, "stddev")
else
n select= casesize(cmdpopulation)
call casemove_alloc("stddev"population, tmp)
stddev = sqrt(stat_object(0.0, "variance"))
case("variance")
m = stat_object(0.0, "mean")
if ( num > 0 ) then
stddev = summa2/real(num) - m*m
else
stddev = 0.0
end if
case("count")
stddev = real(num)
case("mean")
if ( num > 0 ) then
stddev = summa/real(num)
else
stddev = 0.0
end if
case("reset")
summa = 0.0
summa2 = 0.0
num = 0
case default
stddev = 0.0
end select
end if
 
end function stat_object
 
end program Test_Stddev</lang>
 
===Using built-in array awareness===
{{incorrect|Fortran|Function takes array of data at once instead of single values at a time.}}
This uses Fortran's built-in array features (which aren't available in C)
 
{{works with|Fortran|95 and later}}
 
<lang fortran>
program stats
implicit none
 
integer, parameter :: N = 8
integer :: data(N)
real(8) :: mean
real(8) :: std_dev1, std_dev2
 
! Set the data
data = [2,4,4,4,5,5,7,9] ! Strictly this is a Fortran 2003 construct
 
! Use intrinsic function 'sum' to calculate the mean
mean = sum(data)/N
 
! Method1:
! Calculate the standard deviation directly from the definition
std_dev1 = sqrt(sum((data - mean)**2)/N)
 
allocate(population(n + 1))
! Method 2:
population(1:n) = tmp
! Use the alternative version that is less susceptible to rounding error
population(n + 1) = val
std_dev2 = sqrt(sum(data**2)/N - mean**2)
endif
end subroutine sample_add
 
! Calculates standard deviation for given set of values
write(*,'(a,8i2)') 'Data = ',data
real(kind=dp) function stddev(vals)
write(*,'(a,f3.1)') 'Mean = ',mean
real(kind=dp), dimension(:), intent(in) :: vals
write(*,'(a,f3.1)') 'Standard deviation (method 1) = ',std_dev1
real(kind=dp) :: mean
write(*,'(a,f3.1)') 'Standard deviation (method 2) = ',std_dev2
integer(kind=4) :: n
 
n = size(vals)
end program stats
mean = sum(vals)/n
stddev = sqrt(sum((vals - mean)**2)/n)
end function stddev
end program standard_deviation
</lang>
 
Anonymous user
Cookies help us deliver our services. By using our services, you agree to our use of cookies.