Cumulative standard deviation
You are encouraged to solve this task according to the task description, using any language you may know.
Write a stateful function, class, generator or coroutine that takes a series of floating point numbers, one at a time, and returns the running standard deviation of the series. The task implementation should use the most natural programming style of those listed for the function in the implementation language; the task must state which is being used. Do not apply Bessel's Correction; the returned standard deviation should always be computed as if the sample seen so far is the entire population.
Use this to compute the standard deviation of the demonstration set:
C
Of course this function is not thread-safe nor it can be used to compute the standard deviation just for a set of values per time.
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
enum Action { VALUE, STDDEV, MEAN, VAR, COUNT, RESET }; double stat_object(double v, enum Action action) {
static double sum = 0.0; static double sum2 = 0.0; static size_t num = 0;
double m;
switch(action) { case VALUE: num++; sum += v; sum2 += v*v; return stat_object(0.0, STDDEV); case STDDEV: return sqrt(stat_object(0.0, VAR)); case MEAN: return (num>0) ? sum/(double)num : 0.0; case VAR: m = stat_object(0.0, MEAN); return (num>0) ? (sum2/(double)num - m*m) : 0.0; case COUNT: return num; case RESET: sum = sum2 = 0.0; num = 0; return 0.0; }
}</lang>
<lang c>double v[] = { 2,4,4,4,5,5,7,9 };
int main() {
int i; double sd;
for(i=0; i < sizeof(v)/sizeof(double) ; i++) sd = stat_object(v[i], VALUE);
printf("standard dev = %lf\n", sd);
return 0;
}</lang>
Forth
<lang forth>
- st-count ( stats -- n ) f@ ;
- st-sum ( stats -- sum ) float+ f@ ;
- st-sumsq ( stats -- sum*sum ) 2 floats + f@ ;
- st-add ( fnum stats -- )
dup f@ 1e f+ dup f! float+ dup f@ fover f+ dup f! float+ fdup f* dup f@ f+ f! ;
- st-mean ( stats -- mean )
dup st-sum st-count f/ ;
- st-variance ( stats -- var )
dup st-sumsq dup st-mean fdup f* dup st-count f* f- st-count f/ ;
- st-stddev ( stats -- stddev )
st-variance fsqrt ;
create stats 0e f, 0e f, 0e f,
2e stats st-add 4e stats st-add 4e stats st-add 4e stats st-add 5e stats st-add 5e stats st-add 7e stats st-add 9e stats st-add
stats st-stddev f. \ 2. </lang>
Fortran
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.
<lang fortran>program Test_Stddev
implicit none
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
print *, "std dev = ", sd
contains
recursive function stat_object(a, cmd) result(stddev) real :: stddev real, intent(in) :: a character(len=*), intent(in), optional :: cmd
real, save :: summa = 0.0, summa2 = 0.0 integer, save :: num = 0
real :: m
if ( .not. present(cmd) ) then num = num + 1 summa = summa + a summa2 = summa2 + a*a stddev = stat_object(0.0, "stddev") else select case(cmd) case("stddev") 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>
Objective-C
<lang objc>#import <stdio.h>
- import <math.h>
- import <objc/Object.h>
@interface SDAccum : Object {
double sum, sum2; unsigned int num;
} -(id)init; -(double)value: (double)v; -(unsigned int)count; -(double)mean; -(double)variance; -(double)stddev; @end
@implementation SDAccum -(id)init {
sum = sum2 = 0.0; num = 0; return self;
} -(double)value: (double)v {
sum = sum + v; sum2 = sum2 + v*v; num++; return [self stddev];
} -(unsigned int)count {
return num;
} -(double)mean {
return (num>0) ? sum/(double)num : 0.0;
} -(double)variance {
double m = [self mean]; return (num>0) ? (sum2/(double)num - m*m) : 0.0;
} -(double)stddev {
return sqrt([self variance]);
} @end</lang>
<lang objc>double v[] = { 2,4,4,4,5,5,7,9 };
int main() {
int i; double sd;
id sdacc = [SDAccum new];
for(i=0; i < sizeof(v)/sizeof(double) ; i++) sd = [sdacc value: v[i]];
printf("std dev = %lf\n", sd);
return 0;
}</lang>
Perl
<lang perl>{
package SDAccum; sub new {
my $class = shift; my $self = {}; $self->{sum} = 0.0; $self->{sum2} = 0.0; $self->{num} = 0; bless $self, $class; return $self;
} sub count {
my $self = shift; return $self->{num};
} sub mean {
my $self = shift; return ($self->{num}>0) ? $self->{sum}/$self->{num} : 0.0;
} sub variance {
my $self = shift; my $m = $self->mean; return ($self->{num}>0) ? $self->{sum2}/$self->{num} - $m * $m : 0.0;
} sub stddev {
my $self = shift; return sqrt($self->variance);
} sub value {
my $self = shift; my $v = shift; $self->{sum} += $v; $self->{sum2} += $v * $v; $self->{num}++; return $self->stddev;
}
}</lang>
<lang perl>my $sdacc = SDAccum->new; my $sd;
foreach my $v ( 2,4,4,4,5,5,7,9 ) {
$sd = $sdacc->value($v);
} print "std dev = $sd\n";</lang>
Smalltalk
<lang smalltalk>Object subclass: SDAccum [
|sum sum2 num| SDAccum class >> new [ |o| o := super basicNew. ^ o init. ] init [ sum := 0. sum2 := 0. num := 0 ] value: aValue [ sum := sum + aValue. sum2 := sum2 + ( aValue * aValue ). num := num + 1. ^ self stddev ] count [ ^ num ] mean [ num>0 ifTrue: [^ sum / num] ifFalse: [ ^ 0.0 ] ] variance [ |m| m := self mean. num>0 ifTrue: [^ (sum2/num) - (m*m) ] ifFalse: [ ^ 0.0 ] ] stddev [ ^ (self variance) sqrt ]
].</lang>
<lang smalltalk>|sdacc sd| sdacc := SDAccum new.
- ( 2 4 4 4 5 5 7 9 ) do: [ :v | sd := sdacc value: v ].
('std dev = %1' % { sd }) displayNl.</lang>
Tcl
With a Class
<lang tcl>oo::class create SDAccum {
variable sum sum2 num constructor {} { set sum 0.0 set sum2 0.0 set num 0 } method value {x} { set sum2 [expr {$sum2 + $x**2}] set sum [expr {$sum + $x}] incr num return [my stddev] } method count {} { return $num } method mean {} { expr {$sum / double($num)} } method variance {} { expr {($sum2 - double($num)*[my mean]**2)/double($num)} } method stddev {} { expr {sqrt([my variance])} }
}
- Demonstration
set sdacc [SDAccum new] foreach val {2 4 4 4 5 5 7 9} {
set sd [$sdacc value $val]
} puts "the standard deviation is: $sd"</lang> which produces the output:
the standard deviation is: 2.0