Cumulative standard deviation

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 this demonstration set, ${\displaystyle \{2,4,4,4,5,5,7,9\}}$, which is ${\displaystyle 2}$.

procedure Test_Deviation is

```  type Sample is record
N       : Natural := 0;
Mean    : Float := 0.0;
Squares : Float := 0.0;
end record;
procedure Add (Data : in out Sample; Point : Float) is
begin
Data.N       := Data.N + 1;
Data.Mean    := Data.Mean    + Point;
Data.Squares := Data.Squares + Point ** 2;
function Deviation (Data : Sample) return Float is
begin
return Sqrt (Data.Squares / Float (Data.N) - (Data.Mean / Float (Data.N)) ** 2);
end Deviation;
```
```  Data : Sample;
Test : array (1..8) of Float := (2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0);
```

begin

```  for Item in Test'Range loop
end loop;
Put_Line ("Deviation" & Float'Image (Deviation (Data)));
```

end Test_Deviation;</lang> Sample output:

`Deviation 2.00000E+00`

ALGOL 68

Translation of: C
Works with: ALGOL 68 version Standard - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

Note: the use of a UNION to mimic C's enumerated types is "experimental" and probably not typical of "production code". However it is a example of ALGOL 68s conformity CASE clause useful for classroom dissection. <lang Algol68>MODE VALUE = STRUCT(CHAR value),

```    STDDEV = STRUCT(CHAR stddev),
MEAN = STRUCT(CHAR mean),
VAR = STRUCT(CHAR var),
COUNT = STRUCT(CHAR count),
RESET = STRUCT(CHAR reset);
```

MODE ACTION = UNION ( VALUE, STDDEV, MEAN, VAR, COUNT, RESET );

LONG REAL sum := 0; LONG REAL sum2 := 0; INT num := 0;

PROC stat object = (LONG REAL v, ACTION action)LONG REAL: (

``` LONG REAL m;

CASE action IN
(VALUE):(
num +:= 1;
sum +:= v;
sum2 +:= v*v;
stat object(0, LOC STDDEV)
),
(STDDEV):
long sqrt(stat object(0, LOC VAR)),
(MEAN):
IF num>0 THEN sum/LONG REAL(num) ELSE 0 FI,
(VAR):(
m := stat object(0, LOC MEAN);
IF num>0 THEN sum2/LONG REAL(num)-m*m ELSE 0 FI
),
(COUNT):
num,
(RESET):
sum := sum2 := num := 0
ESAC
```

);

[]LONG REAL v = ( 2,4,4,4,5,5,7,9 );

main: (

``` LONG REAL sd;

FOR i FROM LWB v TO UPB v DO
sd := stat object(v[i], LOC VALUE)
OD;

printf((\$"standard dev := "g(0,6)l\$, sd))
```

)</lang>Output:

```standard dev := 2.000000
```
Translation of: python
Works with: ALGOL 68 version Standard - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

A code sample in an object oriented style: <lang Algol68>MODE STAT = STRUCT(

``` LONG REAL sum,
LONG REAL sum2,
INT num
```

);

OP INIT = (REF STAT new)REF STAT:

``` (init OF class stat)(new);
```

MODE CLASSSTAT = STRUCT(

``` PROC (REF STAT, LONG REAL #value#)VOID plusab,
PROC (REF STAT)LONG REAL stddev, mean, variance, count,
PROC (REF STAT)REF STAT init
```

);

CLASSSTAT class stat;

plusab OF class stat := (REF STAT self, LONG REAL value)VOID:(

```   num OF self +:= 1;
sum OF self +:= value;
sum2 OF self +:= value*value
);
```

OP +:= = (REF STAT lhs, LONG REAL rhs)VOID: # some syntatic sugar #

``` (plusab OF class stat)(lhs, rhs);
```

stddev OF class stat := (REF STAT self)LONG REAL:

```   long sqrt((variance OF class stat)(self));
```

OP STDDEV = ([]LONG REAL value)LONG REAL: ( # more syntatic sugar #

``` REF STAT stat = INIT LOC STAT;
FOR i FROM LWB value TO UPB value DO
stat +:= value[i]
OD;
(stddev OF class stat)(stat)
```

);

mean OF class stat := (REF STAT self)LONG REAL:

```   sum OF self/LONG REAL(num OF self);
```

variance OF class stat := (REF STAT self)LONG REAL:(

```   LONG REAL m = (mean OF class stat)(self);
sum2 OF self/LONG REAL(num OF self)-m*m
);
```

count OF class stat := (REF STAT self)LONG REAL:

```   num OF self;
```

init OF class stat := (REF STAT self)REF STAT:(

```   sum OF self := sum2 OF self := num OF self := 0;
self
);
```

[]LONG REAL value = ( 2,4,4,4,5,5,7,9 );

main: (

``` printf((\$"standard deviation operator = "g(0,6)l\$, STDDEV value));
```
``` REF STAT stat = INIT LOC STAT;
```
``` FOR i FROM LWB value TO UPB value DO
stat +:= value[i]
OD;

printf((\$"standard deviation = "g(0,6)l\$, (stddev OF class stat)(stat)));
printf((\$"mean = "g(0,6)l\$, (mean OF class stat)(stat)));
printf((\$"variance = "g(0,6)l\$, (variance OF class stat)(stat)));
printf((\$"count = "g(0,6)l\$, (count OF class stat)(stat)))
```

)</lang>Output:

```standard deviation operator = 2.000000
standard deviation = 2.000000
mean = 5.000000
variance = 4.000000
count = 8.000000
```
Translation of: python
Works with: ALGOL 68 version Standard - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

A simple - but "unpackaged" - code example, useful if the standard deviation is required on only one set of concurrent data: <lang Algol68>LONG REAL sum, sum2; INT n;

PROC sd = (LONG REAL x)LONG REAL:(

```   sum  +:= x;
sum2 +:= x*x;
n    +:= 1;
IF n = 0 THEN 0 ELSE long sqrt(sum2/n - sum*sum/n/n) FI
```

);

sum := sum2 := n := 0; []LONG REAL values = (2,4,4,4,5,5,7,9); FOR i TO UPB values DO

```   LONG REAL value = values[i];
printf((\$2(xg(0,6))l\$, value, sd(value)))
```

OD</lang>Output:

``` 2.000000 .000000
4.000000 1.000000
4.000000 .942809
4.000000 .866025
5.000000 .979796
5.000000 1.000000
7.000000 1.399708
9.000000 2.000000
```

AutoHotkey

ahk forum: discussion <lang AutoHotkey>std(2),std(4),std(4),std(4),std(5),std(5),std(7) MsgBox % std(9) ; 2

std(x="") {

``` Static sum:=0, sqr:=0, n:=0
If (x="")                    ; blank parameter: reset
sum := 0, sqr := 0, n := 0
Else
sum += x, sqr += x*x, n++ ; update state
Return sqrt((sqr-sum*sum/n)/n)
```

}</lang>

Axiom

We implement a domain with dependent type T with the operation + and identity 0:<lang Axiom>)abbrev package TESTD TestDomain TestDomain(T : Join(Field,RadicalCategory)): Exports == Implementation where

``` R ==> Record(n : Integer, sum : T, ssq : T)
Exports == AbelianMonoid with
_+ : (%,T) -> %
_+ : (T,%) -> %
sd : % -> T
Rep := R   -- similar representation and implementation
obj : %
0 == [0,0,0]
obj + (obj2:%) == [obj.n + obj2.n, obj.sum + obj2.sum, obj.ssq + obj2.ssq]
obj + (x:T) == obj + [1, x, x*x]
(x:T) + obj == obj + x
sd obj ==
mean : T := obj.sum / (obj.n::T)
sqrt(obj.ssq / (obj.n::T) - mean*mean)</lang>This can be called using:<lang Axiom>T ==> Expression Integer
```

D ==> TestDomain(T) items := [2,4,4,4,5,5,7,9+x] :: List T; map(sd, scan(+, items, 0\$D))

```                                       +---------------+
+-+  +-+   +-+     +-+  |  2
2\|2  \|3  2\|6    4\|6  \|7x  + 64x + 256
(1)  [0,1,-----,----,-----,1,-----,------------------]
3     2    5       7            8
Type: List(Expression(Integer))
```

eval subst(%.8,x=0)

```   (2)  2
Type: Expression(Integer)</lang>
```

BBC BASIC

Uses the MOD(array()) and SUM(array()) functions. <lang bbcbasic> MAXITEMS = 100

```     FOR i% = 1 TO 8
PRINT "Value = "; n ", running SD = " FNrunningsd(n)
NEXT
END

DATA 2,4,4,4,5,5,7,9

DEF FNrunningsd(n)
PRIVATE list(), i%
DIM list(MAXITEMS)
i% += 1
list(i%) = n
= SQR(MOD(list())^2/i% - (SUM(list())/i%)^2)</lang>
```

Output:

```Value = 2, running SD = 0
Value = 4, running SD = 1
Value = 4, running SD = 0.942809043
Value = 4, running SD = 0.866025404
Value = 5, running SD = 0.979795901
Value = 5, running SD = 1
Value = 7, running SD = 1.39970842
Value = 9, running SD = 2
```

C

<lang c>#include <stdio.h>

1. include <stdlib.h>
2. include <math.h>

enum Action { STDDEV, MEAN, VAR, COUNT };

typedef struct stat_obj_struct {

```  double sum, sum2;
size_t num;
Action action;
```

} sStatObject, *StatObject;

StatObject NewStatObject( Action action ) {

``` so = malloc(sizeof(sStatObject));
so->sum = 0.0;
so->sum2 = 0.0;
so->num = 0;
so->action = action;
return so;
```

}

1. define FREE_STAT_OBJECT(so) \
```  free(so); so = NULL
```

double stat_obj_value(StatObject so, Action action) {

``` double num, mean, var, stddev;

if (so->num == 0.0) return 0.0;
num = so->num;
if (action==COUNT) return num;
mean = so->sum/num;
if (action==MEAN) return mean;
var = so->sum2/num - mean*mean;
if (action==VAR) return var;
stddev = sqrt(var);
if (action==STDDEV) return stddev;
return 0;
```

}

double stat_object_add(StatObject so, double v) {

``` so->num++;
so->sum += v;
so->sum2 += v*v;
return stat_obj_value(so->action);
```

}</lang>

<lang c>double v[] = { 2,4,4,4,5,5,7,9 };

int main() {

``` int i;
StatObject so = NewStatObject( STDDEV );
```
``` for(i=0; i < sizeof(v)/sizeof(double) ; i++)
printf("val: %lf  std dev: %lf\n", v[i], stat_object_add(so, v[i]));
```
``` FREE_STAT_OBJECT(so);
return 0;
```

}</lang>

C++

<lang cpp>#include <algorithm>

1. include <iostream>
2. include <iterator>
3. include <cmath>
4. include <vector>
5. include <iterator>
6. include <numeric>

template <typename Iterator> double standard_dev( Iterator begin , Iterator end ) {

```  double mean = std::accumulate( begin , end , 0 ) / std::distance( begin , end ) ;
std::vector<double> squares ;
for( Iterator vdi = begin ; vdi != end ; vdi++ )
squares.push_back( std::pow( *vdi - mean , 2 ) ) ;
return std::sqrt( std::accumulate( squares.begin( ) , squares.end( ) , 0 ) / squares.size( ) ) ;
```

}

int main( ) {

```  double demoset[] = { 2 , 4 , 4 , 4 , 5 , 5 , 7 , 9 } ;
int demosize = sizeof demoset / sizeof *demoset ;
std::cout << "The standard deviation of\n" ;
std::copy( demoset , demoset + demosize , std::ostream_iterator<double>( std::cout, " " ) ) ;
std::cout << "\nis " << standard_dev( demoset , demoset + demosize ) << " !\n" ;
return 0 ;
```

}</lang>

C#

<lang csharp>using System; using System.Collections.Generic; using System.Linq;

namespace standardDeviation {

```   class Program
{
static void Main(string[] args)
{
List<double> nums = new List<double> { 2, 4, 4, 4, 5, 5, 7, 9 };
for (int i = 1; i <= nums.Count; i++)
Console.WriteLine(sdev(nums.GetRange(0, i)));
}
```
```       static double sdev(List<double> nums)
{
List<double> store = new List<double>();
foreach (double n in nums)
store.Add((n - nums.Average()) * (n - nums.Average()));
```
```           return Math.Sqrt(store.Sum() / store.Count);
}
}
```

}</lang>

```0
1
0,942809041582063
0,866025403784439
0,979795897113271
1
1,39970842444753
2```

Clojure

<lang lisp> (defn std-dev [samples]

``` (let [n (count samples)
```

mean (/ (reduce + samples) n) intermediate (map #(Math/pow (- %1 mean) 2) samples)]

```   (Math/sqrt
(/ (reduce + intermediate) n))))
```

(println (std-dev [2 4 4 4 5 5 7 9])) ;;2.0

</lang>

Common Lisp

<lang lisp>(defun std-dev (samples)

``` (let* ((n (length samples))
```

(mean (/ (reduce #'+ samples) n)) (tmp (mapcar (lambda (x) (sqr (- x mean))) samples)))

```   (sqrt (/ (reduce #'+ tmp) n))))
```

(format t "~a" (std-dev '(2 4 4 4 5 5 7 9))) </lang>

Based on some googled web site; written ages ago.

<lang lisp>(defun arithmetic-average (samples)

``` (/ (reduce #'+ samples)
(length samples)))
```

(defun standard-deviation (samples)

``` (let ((mean (arithmetic-average samples)))
(sqrt (* (/ 1.0d0 (length samples))
(reduce #'+ samples
:key (lambda (x)
(expt (- x mean) 2)))))))
```

(defun make-deviator ()

``` (let ((numbers '()))
(lambda (x)
(push x numbers)
(standard-deviation numbers))))</lang>
```

<lang lisp>CL-USER> (loop with deviator = (make-deviator)

```              for i in '(2 4 4 4 5 5 7 9)
collect (list i (funcall deviator i))) ==>
```

((2 0.0d0)

```(4 1.0d0)
(4 0.9428090415820634d0)
(4 0.8660254037844386d0)
(5 0.9797958971132713d0)
(5 1.0d0)
(7 1.3997084244475304d0)
(9 2.0d0))</lang>
```

Since we don't care about the sample values once std dev is computed, we only need to keep track of their sum and square sums, hence:<lang lisp>(defun running-stddev ()

``` (let ((sum 0) (sq 0) (n 0))
(lambda (x)
(incf sum x) (incf sq (* x x)) (incf n)
(/ (sqrt (- (* n sq) (* sum sum))) n))))
```

(loop with f = (running-stddev) for i in '(2 4 4 4 5 5 7 9) do (format t "~a ~a~%" i (funcall f i)))</lang>

D

<lang d>import std.stdio, std.math;

struct StdDev {

```   real sum = 0.0, sqSum = 0.0;
long nvalues;
```
```   void addNumber(in real input) pure nothrow {
nvalues++;
sum += input;
sqSum += input ^^ 2;
}
```
```   real getStdDev() const pure nothrow {
if (nvalues == 0)
return 0.0;
immutable real mean = sum / nvalues;
return sqrt(sqSum / nvalues - mean ^^ 2);
}
```

}

void main() {

```   StdDev stdev;
```
```   foreach (el; [2.0, 4, 4, 4, 5, 5, 7, 9]) {
writefln("%e", stdev.getStdDev());
}
```

}</lang> Output:

```0.000000e+00
1.000000e+00
9.428090e-01
8.660254e-01
9.797959e-01
1.000000e+00
1.399708e+00
2.000000e+00```

Delphi

Delphi has 2 functions in Math unit for standard deviation: StdDev (with Bessel correction) and PopnStdDev (without Bessel correction). The task assumes the second function:

<lang Delphi>program StdDevTest;

{\$APPTYPE CONSOLE}

uses

``` Math;
```

begin

``` Writeln(PopnStdDev([2,4,4,4,5,5,7,9]));
```

end.</lang>

E

This implementation produces two (function) objects sharing state. It is idiomatic in E to separate input from output (read from write) rather than combining them into one object.

The algorithm is

Translation of: Perl

and the results were checked against #Python.

<lang e>def makeRunningStdDev() {

```   var sum := 0.0
var sumSquares := 0.0
var count := 0.0

def insert(v) {
sum += v
sumSquares += v ** 2
count += 1
}

/** Returns the standard deviation of the inputs so far, or null if there
have been no inputs. */
def stddev() {
if (count > 0) {
def meanSquares := sumSquares/count
def mean := sum/count
def variance := meanSquares - mean**2
return variance.sqrt()
}
}

return [insert, stddev]
```

}</lang>

<lang e>? def [insert, stddev] := makeRunningStdDev()

1. value: <insert>, <stddev>

? [stddev()]

1. value: [null]

? for value in [2,4,4,4,5,5,7,9] { > insert(value) > println(stddev()) > } 0.0 1.0 0.9428090415820626 0.8660254037844386 0.9797958971132716 1.0 1.3997084244475297 2.0</lang>

F#

<lang fsharp>module StdDev =

```   let main =
let running_stddev (N, M, M2, SD) x =
let n = N + 1.
let delta = x - M
let mean = M + delta/n
let m2 = M2 + delta * (x - mean)
let sd = (sqrt (m2/n))
printfn "%.1f mean:%.4f sd:%f" x mean sd
(n, mean, m2, sd)
let (count, mean, state, sd) =
[ 2.; 4.; 4.; 4.; 5.; 5.; 7.; 9. ]
|> List.fold (fun (a,b,c,d) v -> running_stddev (a,b,c,d) v) (0., 0., 0., 0.)
printfn "stddev: %f" sd
```

do StdDev.main</lang>

Output:

```2.0 mean:2.0000 sd:0.000000
4.0 mean:3.0000 sd:1.000000
4.0 mean:3.3333 sd:0.942809
4.0 mean:3.5000 sd:0.866025
5.0 mean:3.8000 sd:0.979796
5.0 mean:4.0000 sd:1.000000
7.0 mean:4.4286 sd:1.399708
9.0 mean:5.0000 sd:2.000000
stddev: 2.000000```

Factor

<lang factor>USING: accessors io kernel math math.functions math.parser sequences ; IN: standard-deviator

TUPLE: standard-deviator sum sum^2 n ;

<standard-deviator> ( -- standard-deviator )
```   0.0 0.0 0 standard-deviator boa ;
```
current-std ( standard-deviator -- std )
```   [ [ sum^2>> ] [ n>> ] bi / ]
[ [ sum>> ] [ n>> ] bi / sq ] bi - sqrt ;
```
add-value ( value standard-deviator -- )
```   [ nip [ 1 + ] change-n drop ]
[ [ + ] change-sum drop ]
[ [ [ sq ] dip + ] change-sum^2 drop ] 2tri ;
```
main ( -- )
```   { 2 4 4 4 5 5 7 9 }
<standard-deviator> [ [ add-value ] curry each ] keep
current-std number>string print ;</lang>
```

Forth

<lang forth>: f+! ( x addr -- ) dup f@ f+ f! ;

st-count ( stats -- n ) f@ ;
st-sum ( stats -- sum ) float+ f@ ;
st-sumsq ( stats -- sum*sum ) 2 floats + f@ ;
st-add ( fnum stats -- )
```   1e dup f+!  float+
fdup dup f+!  float+
fdup 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 ;</lang>
```

This variation is more numerically stable when there are large numbers of samples or large sample ranges. <lang forth>: st-count ( stats -- n ) f@ ;

st-mean ( stats -- mean ) float+ f@ ;
st-nvar ( stats -- n*var ) 2 floats + f@ ;
st-variance ( stats -- var ) dup st-nvar st-count f/ ;
st-stddev ( stats -- stddev ) st-variance fsqrt ;
st-add ( x stats -- )
``` 1e dup f+!			\ update count
fdup dup st-mean f- fswap
( delta x )
fover dup st-count f/
( delta x delta/n )
float+ dup f+!		\ update mean
( delta x )
dup f@ f-  f*  float+ f+! ;	\ update nvar</lang>
```

<lang forth>create stats 0e f, 0e f, 0e f,

stats st-stddev f. \ 2.</lang>

Fortran

Translation of: C
Works with: Fortran version 95 and later

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>

Using built-in array awareness

This uses Fortran's built-in array features (which aren't available in C)

Works with: Fortran version 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)
```
``` ! Method 2:
! Use the alternative version that is less susceptible to rounding error
std_dev2 = sqrt(sum(data**2)/N - mean**2)
```
``` write(*,'(a,8i2)') 'Data = ',data
write(*,'(a,f3.1)') 'Mean = ',mean
write(*,'(a,f3.1)') 'Standard deviation (method 1) = ',std_dev1
write(*,'(a,f3.1)') 'Standard deviation (method 2) = ',std_dev2
```

end program stats </lang>

Go

Algorithm to reduce rounding errors from WP article.

State maintained with a closure. <lang go>package main

import (

```   "fmt"
"math"
```

)

func newRsdv() func(float64) float64 {

```   var n, a, q  float64
return func(x float64) float64 {
n++
a1 := a+(x-a)/n
q, a = q+(x-a)*(x-a1), a1
return math.Sqrt(q/n)
}
```

}

func main() {

```   r := newRsdv()
for _, x := range []float64{2,4,4,4,5,5,7,9} {
fmt.Println(r(x))
}
```

}</lang> Output:

```0
1
0.9428090415820634
0.8660254037844386
0.9797958971132713
1
1.3997084244475304
2
```

Groovy

Solution: <lang groovy>def sum = 0 def sumSq = 0 def count = 0 [2,4,4,4,5,5,7,9].each {

```   sum += it
sumSq += it*it
count++
println "running std.dev.: \${(sumSq/count - (sum/count)**2)**0.5}"
```

}</lang>

Output:

```running std.dev.: 0
running std.dev.: 1
running std.dev.: 0.9428090416999145
running std.dev.: 0.8660254037844386
running std.dev.: 0.9797958971132712
running std.dev.: 1
running std.dev.: 1.3997084243469262
running std.dev.: 2```

We store the state in the `ST` monad using an `STRef`.

sd :: [Double] -> Double sd l = sqrt \$ sum (map ((^2) . subtract mean) l) / n

``` where n = genericLength l
mean = sum l / n
```

sdAccum :: ST s (Double -> ST s Double) sdAccum = do

```   accum <- newSTRef []
return \$ \x -> do
modifySTRef accum (x:)
return \$ sd list
```

main = mapM_ print results

``` where results = runST \$ do
runningSD <- sdAccum
mapM runningSD [2, 4, 4, 4, 5, 5, 7, 9]</lang>
```

Haxe

<lang haxe>using Lambda;

class Main { static function main():Void { var nums = [2, 4, 4, 4, 5, 5, 7, 9]; for (i in 1...nums.length+1) Sys.println(sdev(nums.slice(0, i))); }

static function average<T:Float>(nums:Array<T>):Float { return nums.fold(function(n, t) return n + t, 0) / nums.length; }

static function sdev<T:Float>(nums:Array<T>):Float { var store = []; var avg = average(nums); for (n in nums) { store.push((n - avg) * (n - avg)); }

return Math.sqrt(average(store)); } }</lang>

```0
1
0.942809041582063
0.866025403784439
0.979795897113271
1
1.39970842444753
2```

HicEst

<lang HicEst>REAL :: n=8, set(n), sum=0, sum2=0

set = (2,4,4,4,5,5,7,9)

DO k = 1, n

```  WRITE() 'Adding ' // set(k) // 'stdev = ' // stdev(set(k))
```

ENDDO

END ! end of "main"

FUNCTION stdev(x)

```  USE : sum, sum2, k
sum = sum + x
sum2 = sum2 + x*x
stdev = ( sum2/k - (sum/k)^2) ^ 0.5
END</lang>
```
```Adding 2 stdev = 0

Icon and Unicon

<lang Icon>rocedure main()

stddev() # reset state / empty every s := stddev(![2,4,4,4,5,5,7,9]) do

```  write("stddev (so far) := ",s)
```

end

procedure stddev(x) /: running standard deviation static X,sumX,sum2X

```  if /x then {   # reset state
X := []
sumX := sum2X := 0.
}
else {         # accumulate
put(X,x)
sumX +:= x
sum2X +:= x^2
return sqrt( (sum2X / *X) - (sumX / *X)^2 )
}
```

end</lang>

Sample output:

```stddev (so far) := 0.0
stddev (so far) := 1.0
stddev (so far) := 0.9428090415820626
stddev (so far) := 0.8660254037844386
stddev (so far) := 0.9797958971132716
stddev (so far) := 1.0
stddev (so far) := 1.39970842444753
stddev (so far) := 2.0```

J

J is block-oriented; it expresses algorithms with the semantics of all the data being available at once. It does not have native lexical closure or coroutine semantics. It is possible to implement these semantics in J; see Moving Average for an example. We will not reprise that here. <lang j> mean=: +/ % #

```  dev=: - mean
stddevP=: [: %:@mean *:@dev          NB. A) 3 equivalent defs for stddevP
stddevP=: [: mean&.:*: dev           NB. B) uses Under (&.:) to apply inverse of *: after mean
stddevP=: %:@(mean@:*: - *:@mean)    NB. C) sqrt of ((mean of squares) - (square of mean))
```

```  stddevP\ 2 4 4 4 5 5 7 9
```

0 1 0.942809 0.866025 0.979796 1 1.39971 2</lang>

Alternatives:
Using verbose names for J primitives. <lang j> of =: @:

```  sqrt   =: %:
sum    =: +/
squares=: *:
data   =: ]
mean   =: sum % #
```
```  stddevP=: sqrt of mean of squares of (data-mean)
```
```  stddevP\ 2 4 4 4 5 5 7 9
```

0 1 0.942809 0.866025 0.979796 1 1.39971 2</lang>

Translation of: R

Or we could take a cue from the R implementation and reverse the Bessel correction to stddev:

<lang j> require'stats'

```  (%:@:(%~<:)@:# * stddev)\ 2 4 4 4 5 5 7 9
```

0 1 0.942809 0.866025 0.979796 1 1.39971 2</lang>

Java

Works with: Java version 1.5+

<lang java5>import java.util.LinkedList; public class StdDev {

```   LinkedList<Double> nums;
public StdDev() {
}
public static void main(String[] args) {
double[] testData = {2,4,4,4,5,5,7,9};
StdDev sd = new StdDev();

for (double x : testData) {
sd.newNum(x);
System.out.println(sd.getSD());
}
}

public void newNum(double num) {
}

public double getAvg() {
if (nums.isEmpty()) return 0;
double ret = 0;
double sum = 0;
for (double num : nums) {
sum += num;
}
return sum / nums.size();
}

public double getSD() {
double sqDiffs = 0;
double avg = getAvg();
for (double num : nums) {
sqDiffs += (num - avg) * (num - avg);
}
return Math.sqrt(sqDiffs / nums.size());
}
```

}</lang>

JavaScript

Uses a closure. <lang javascript>function running_stddev() {

```   var n = 0;
var sum = 0.0;
var sum_sq = 0.0;
return function(num) {
n++;
sum += num;
sum_sq += num*num;
return Math.sqrt( (sum_sq / n) - Math.pow(sum / n, 2) );
}
```

}

var sd = running_stddev(); var nums = [2,4,4,4,5,5,7,9]; var stddev = []; for (var i in nums)

```   stddev.push( sd(nums[i]) );
```

// using WSH WScript.Echo(stddev.join(', ');</lang>

output:

`0, 1, 0.942809041582063, 0.866025403784439, 0.979795897113273, 1, 1.39970842444753, 2`

Liberty BASIC

Using a global array to maintain the state. Implements definition explicitly. <lang lb>

```   dim SD.storage\$( 100)   '   can call up to 100 versions, using ID to identify.. arrays are global.
'   holds (space-separated) number of data items so far, current sum.of.values and current sum.of.squares
```
```   for i =1 to 8
print "New data "; x; " so S.D. now = "; using( "###.######", standard.deviation( 1, x))
next i
```
```   end
```

function standard.deviation( ID, in)

``` if SD.storage\$( ID) ="" then SD.storage\$( ID) ="0 0 0"
num.so.far =val( word\$( SD.storage\$( ID), 1))
sum.vals   =val( word\$( SD.storage\$( ID), 2))
sum.sqs    =val( word\$( SD.storage\$( ID), 3))
num.so.far =num.so.far +1
sum.vals   =sum.vals   +in
sum.sqs    =sum.sqs    +in^2
```
``` ' standard deviation = square root of (the average of the squares less the square of the average)
standard.deviation   =(               ( sum.sqs /num.so.far)      -    ( sum.vals /num.so.far)^2)^0.5
```
``` SD.storage\$( ID) =str\$( num.so.far) +" " +str\$( sum.vals) +" " +str\$( sum.sqs)
```

end function

```   Data 2, 4, 4, 4, 5, 5, 7, 9
```

</lang>

```New data 2 so S.D. now =   0.000000
New data 4 so S.D. now =   1.000000
New data 4 so S.D. now =   0.942809
New data 4 so S.D. now =   0.866025
New data 5 so S.D. now =   0.979796
New data 5 so S.D. now =   1.000000
New data 7 so S.D. now =   1.399708
New data 9 so S.D. now =   2.000000
```

lua

Uses a closure. Translation of JavaScript. <lang lua>function stdev()

``` local sum, sumsq, k = 0,0,0
return function(n)
sum, sumsq, k = sum + n, sumsq + n^2, k+1
return math.sqrt((sumsq / k) - (sum/k)^2)
end
```

end

ldev = stdev() for i, v in ipairs{2,4,4,4,5,5,7,9} do

``` print(ldev(v))
```

end</lang>

Mathematica

<lang Mathematica>runningSTDDev[n_] := (If[Not[ValueQ[\$Data]], \$Data = {}];

``` StandardDeviation[AppendTo[\$Data, n]])</lang>
```

MATLAB / Octave

The simple form is, computing only the standand deviation of the whole data set:

<lang Matlab> x = [2,4,4,4,5,5,7,9];

``` n = length (x);
```
``` m  = mean (x);
x2 = mean (x .* x);
dev= sqrt (x2 - m * m)
dev = 2 </lang>
```

When the intermediate results are also needed, one can use this vectorized form:

<lang Matlab> m = cumsum(x) ./ [1:n]; % running mean

``` x2= cumsum(x.^2) ./ [1:n];   % running squares
```
``` dev = sqrt(x2 - m .* m)
dev =
0.00000   1.00000   0.94281   0.86603   0.97980   1.00000   1.39971   2.00000
```

</lang>

Objective-C

<lang objc>#import <Foundation/Foundation.h>

@interface SDAccum : NSObject {

``` double sum, sum2;
unsigned int num;
```

} -(double)value: (double)v; -(unsigned int)count; -(double)mean; -(double)variance; -(double)stddev; @end

@implementation SDAccum -(double)value: (double)v {

``` sum += v;
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

int main() {

``` double v[] = { 2,4,4,4,5,5,7,9 };
```
``` SDAccum *sdacc = [SDAccum new];
```
``` for(int i=0; i < sizeof(v)/sizeof(*v) ; i++)
printf("adding %f\tstddev = %f\n", v[i], [sdacc value: v[i]]);
```
``` [sdacc release];
```
``` return 0;
```

}</lang>

Blocks

Works with: Mac OS X version 10.6+
Works with: iOS version 4+

<lang objc>#import <Foundation/Foundation.h>

typedef double (^Func)(double); // a block that takes a double and returns a double

Func sdCreator() {

``` __block int n = 0;
__block double sum = 0;
__block double sum2 = 0;
return [[^(double x) {
sum += x;
sum2 += x*x;
n++;
return sqrt(sum2/n - sum*sum/n/n);
} copy] autorelease];
```

}

int main() {

``` NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
```
``` double v[] = { 2,4,4,4,5,5,7,9 };
```
``` Func sdacc = sdCreator();
```
``` for(int i=0; i < sizeof(v)/sizeof(*v) ; i++)
printf("adding %f\tstddev = %f\n", v[i], sdacc(v[i]));
```
``` [pool release];
return 0;
```

}</lang>

Objeck

Translation of: Java

<lang objeck> use Structure;

bundle Default {

``` class StdDev {
nums : FloatVector;

New() {
nums := FloatVector->New();
}

function : Main(args : String[]) ~ Nil {
sd := StdDev->New();
test_data := [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
each(i : test_data) {
sd->GetSD()->PrintLine();
};
}

method : public : AddNum(num : Float) ~ Nil {
}

method : public : native : GetSD() ~ Float {
sq_diffs := 0.0;
avg := nums->Average();
each(i : nums) {
num := nums->Get(i);
sq_diffs += (num - avg) * (num - avg);
};

return (sq_diffs / nums->Size())->SquareRoot();
}
}
```

} </lang>

OCaml

<lang ocaml>let sqr x = x *. x

let stddev l =

``` let n, sx, sx2 =
List.fold_left
(fun (n, sx, sx2) x -> succ n, sx +. x, sx2 +. sqr x)
(0, 0., 0.) l
in
sqrt ((sx2 -. sqr sx /. float n) /. float n)
```

let _ =

``` let l = [ 2.;4.;4.;4.;5.;5.;7.;9. ] in
Printf.printf "List: ";
List.iter (Printf.printf "%g  ") l;
Printf.printf "\nStandard deviation: %g\n" (stddev l)</lang>
```

Sample output:

```List: 2  4  4  4  5  5  7  9
Standard deviation: 2
```

PARI/GP

Uses the Cramer-Young updating algorithm. For demonstration it displays the mean and variance at each step. <lang parigp>newpoint(x)={

``` myT=x;
myS=0;
myN=1;
[myT,myS]/myN
```

``` myT+=x;
myN++;
myS+=(myN*x-myT)^2/myN/(myN-1);
[myT,myS]/myN
```

``` print(newpoint(v[1]));
print("Mean: ",myT/myN);
print("Standard deviation: ",sqrt(myS/myN))
```

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>

A much shorter version using a closure and a property of the variance:

<lang perl># <(x - <x>)²> = <x²> - <x>² sub stddev; {

```   my @sum;
sub stddev {
```

my \$x = shift; \$sum[0]++; return sqrt( (\$sum[2] += \$x**2) / \$sum[0] - ((\$sum[1] += \$x) / \$sum[0])**2 );

```   }
```

}

print stddev(\$_), "\n" for qw(2 4 4 4 5 5 7 9);</lang>

Output:
```0
1
0.942809041582063
0.866025403784439
0.979795897113272
1
1.39970842444753
2```

Perl 6

Works with: Rakudo Star version 2010.08

Using a closure: <lang perl6>sub sd (@a) {

```   my \$mean = ([+] @a) / @a;
sqrt ([+] map { (\$^x - \$mean)**2 }, @a) / @a;
```

}

sub sdaccum {

```   my @a;
return { push @a, \$^x; sd @a; };
```

};

my &f = sdaccum; say f \$_ for 2, 4, 4, 4, 5, 5, 7, 9;</lang>

Using a state variable: <lang perl6># remember that <(x-<x>)²> = <x²> - <x>² sub stddev(\$x) {

```   state @sum;
return sqrt(
(@sum[2] += \$x**2) / ++@sum[0]   -
((@sum[1] += \$x) / @sum[0])**2
);
```

}

say stddev \$_ for <2 4 4 4 5 5 7 9>;</lang>

Output:
```0
1
0.942809041582063
0.866025403784439
0.979795897113271
1
1.39970842444753
2```

PL/I

<lang PL/I> declare A(200) float; declare ap fixed binary; declare i fixed binary;

ap = 0; do i = 1 to 10;

```  ap = ap + 1;
get list ( A(ap) );
put skip list ( std_dev (a, ap) );
```

end;

std_dev: procedure (A, ap) returns (float);

```  declare A(*) float, ap fixed binary nonassignable;
declare i fixed binary;
declare B(ap) float, average float;
```
```  do i = 1 to ap;
B(i) = A(i);
end;
average = sum(A) /ap;
```
```  return ( sqrt (sum(B**2)/ap - average**2) );
```

end std_dev; </lang>

PicoLisp

<lang PicoLisp>(scl 2)

(de stdDev ()

```  (curry ((Data)) (N)
(push 'Data N)
(let (Len (length Data)  M (*/ (apply + Data) Len))
(sqrt
(*/
(sum
'((N) (*/ (- N M) (- N M) 1.0))
Data )
1.0
Len )
T ) ) ) )
```

(let Fun (stdDev)

```  (for N (2.0 4.0 4.0 4.0 5.0 5.0 7.0 9.0)
(prinl (format N *Scl) " -> " (format (Fun N) *Scl)) ) )</lang>
```

Output:

```2.00 -> 0.00
4.00 -> 1.00
4.00 -> 0.94
4.00 -> 0.87
5.00 -> 0.98
5.00 -> 1.00
7.00 -> 1.40
9.00 -> 2.00```

PowerShell

This implementation takes the form of an advanced function which can act like a cmdlet and receive input from the pipeline. <lang powershell>function Get-StandardDeviation {

```   begin {
\$avg = 0
\$nums = @()
}
process {
\$nums += \$_
\$avg = (\$nums | Measure-Object -Average).Average
\$sum = 0;
\$nums | ForEach-Object { \$sum += (\$avg - \$_) * (\$avg - \$_) }
[Math]::Sqrt(\$sum / \$nums.Length)
}
```

}</lang> Usage as follows:

```PS> 2,4,4,4,5,5,7,9 | Get-StandardDeviation
0
1
0.942809041582063
0.866025403784439
0.979795897113271
1
1.39970842444753
2```

PureBasic

<lang PureBasic>;Define our Standard deviation function Declare.d Standard_deviation(x)

Main program

If OpenConsole()

``` Define i, x
Restore MyList
For i=1 To 8
PrintN(StrD(Standard_deviation(x)))
Next i
Print(#CRLF\$+"Press ENTER to exit"): Input()
```

EndIf

Calculation procedure, with memory

Procedure.d Standard_deviation(In)

``` Static in_summa, antal
in_summa+in
antal+1
```

EndProcedure

data section

DataSection MyList:

``` Data.i  2,4,4,4,5,5,7,9
```

EndDataSection</lang>

Outputs

```0.0000000000
1.0000000000
0.9428090416
0.8660254038
0.9797958971
1.0000000000
1.3997084244
2.0000000000
```

Python

Using a function with attached properties

The program should work with Python 2.x and 3.x, although the output would not be a tuple in 3.x <lang python>>>> from math import sqrt >>> def sd(x):

```   sd.sum  += x
sd.sum2 += x*x
sd.n    += 1.0
sum, sum2, n = sd.sum, sd.sum2, sd.n
return sqrt(sum2/n - sum*sum/n/n)
```

>>> sd.sum = sd.sum2 = sd.n = 0 >>> for value in (2,4,4,4,5,5,7,9):

```   print (value, sd(value))
```

(2, 0.0) (4, 1.0) (4, 0.94280904158206258) (4, 0.8660254037844386) (5, 0.97979589711327075) (5, 1.0) (7, 1.3997084244475311) (9, 2.0) >>></lang>

Using a class instance

<lang python>>>> class SD(object): # Plain () for python 3.x def __init__(self): self.sum, self.sum2, self.n = (0,0,0) def sd(self, x): self.sum += x self.sum2 += x*x self.n += 1.0 sum, sum2, n = self.sum, self.sum2, self.n return sqrt(sum2/n - sum*sum/n/n)

>>> sd_inst = SD() >>> for value in (2,4,4,4,5,5,7,9): print (value, sd_inst.sd(value))</lang>

Using a Closure

Works with: Python version 3.x

<lang python>>>> from math import sqrt >>> def sdcreator(): sum = sum2 = n = 0 def sd(x): nonlocal sum, sum2, n

sum += x sum2 += x*x n += 1.0 return sqrt(sum2/n - sum*sum/n/n) return sd

>>> sd = sdcreator() >>> for value in (2,4,4,4,5,5,7,9): print (value, sd(value))

2 0.0 4 1.0 4 0.942809041582 4 0.866025403784 5 0.979795897113 5 1.0 7 1.39970842445 9 2.0</lang>

R

Built-in Std Dev fn

<lang R>#The built-in standard deviation function applies the Bessel correction. To reverse this, we can apply an uncorrection.

1. If na.rm is true, missing data points (NA values) are removed.
```reverseBesselCorrection <- function(x, na.rm=FALSE)
{
if(na.rm) x <- x[!is.na(x)]
len <- length(x)
if(len < 2) stop("2 or more data points required")
sqrt((len-1)/len)
}
testdata <- c(2,4,4,4,5,5,7,9)
reverseBesselCorrection(testdata)*sd(testdata) #2</lang>
```

From scratch

<lang R>#Again, if na.rm is true, missing data points (NA values) are removed.

```uncorrectedsd <- function(x, na.rm=FALSE)
{
len <- length(x)
if(len < 2) stop("2 or more data points required")
mu <- mean(x, na.rm=na.rm)
ssq <- sum((x - mu)^2, na.rm=na.rm)
usd <- sqrt(ssq/len)
usd
}
uncorrectedsd(testdata) #2</lang>
```

REXX

Works with: oorexx

This is indeed Object REXX.

<lang rexx>sdacc = .SDAccum~new x = .array~of(2,4,4,4,5,5,7,9) sd = 0 do i = 1 to x~size

```  sd = sdacc~value(x[i])
```

end

say "std dev = "sd

class SDAccum
method sum attribute
method sum2 attribute
method count attribute
method init
``` self~sum = 0.0
self~sum2 = 0.0
self~count = 0
```
method value
``` expose sum sum2 count
parse arg x
sum = sum + x
sum2 = sum2 + x*x
count = count + 1
return self~stddev
```
method mean
``` expose sum count
return sum/count
```
method variance
``` expose sum2  count
m = self~mean
return sum2/count - m*m
```
method stddev
``` return self~sqrt(self~variance)
```
method sqrt
``` arg n
if n = 0 then return 0
ans = n / 2
prev = n
do until prev = ans
prev = ans
ans = ( prev + ( n / prev ) ) / 2
end
return ans</lang>
```

Ruby

Object

Uses an object to keep state.

"Simplification of the formula [...] for standard deviation [...] can be memorized as taking the square root of (the average of the squares less the square of the average)." c.f. wikipedia.

<lang ruby>class StdDevAccumulator

``` def initialize
@n, @sum, @sumofsquares = 0, 0.0, 0.0
end

def <<(num)
# return self to make this possible:  sd << 1 << 2 << 3 # => 0.816496580927726
@n += 1
@sum += num
@sumofsquares += num**2
self
end

def stddev
Math.sqrt( (@sumofsquares / @n) - (@sum / @n)**2 )
end

def to_s
stddev.to_s
end
```

end

sd = StdDevAccumulator.new i = 0 [2,4,4,4,5,5,7,9].each {|n| puts "adding #{n}: stddev of #{i+=1} samples is #{sd << n}" }</lang>

```adding 2: stddev of 1 samples is 0.0
adding 4: stddev of 2 samples is 1.0
adding 4: stddev of 3 samples is 0.942809041582063
adding 4: stddev of 4 samples is 0.866025403784439
adding 5: stddev of 5 samples is 0.979795897113272
adding 5: stddev of 6 samples is 1.0
adding 7: stddev of 7 samples is 1.39970842444753
adding 9: stddev of 8 samples is 2.0```

Closure

<lang ruby>def sdaccum

``` n, sum, sum2 = 0, 0.0, 0.0
lambda do |num|
n += 1
sum += num
sum2 += num**2
Math.sqrt( (sum2 / n) - (sum / n)**2 )
end
```

end

sd = sdaccum [2,4,4,4,5,5,7,9].each {|n| print sd.call(n), ", "}</lang>

`0.0, 1.0, 0.942809041582063, 0.866025403784439, 0.979795897113272, 1.0, 1.39970842444753, 2.0, `

Scheme

<lang scheme> (define ((running-stddev . nums) num)

``` (set! nums (cons num nums))
(sqrt (- (/ (apply + (map (lambda (i) (* i i)) nums)) (length nums)) (expt (/ (apply + nums) (length nums)) 2))))
```

</lang>

Smalltalk

Works with: GNU 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.

1. ( 2 4 4 4 5 5 7 9 ) do: [ :v | sd := sdacc value: v ].

('std dev = %1' % { sd }) displayNl.</lang>

Tcl

With a Class

Works with: Tcl version 8.6

or

Library: TclOO

<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 / \$num}
}
method variance {} {
expr {\$sum2/\$num - [my mean]**2}
}
method stddev {} {
expr {sqrt([my variance])}
}
```

}

1. 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`

With a Coroutine

Works with: Tcl version 8.6

<lang tcl># Make a coroutine out of a lambda application coroutine sd apply {{} {

```   set sum 0.0
set sum2 0.0
set sd {}
# Keep processing argument values until told not to...
while {[set val [yield \$sd]] ne "stop"} {
incr n
set sum [expr {\$sum + \$val}]
set sum2 [expr {\$sum2 + \$val**2}]
set sd [expr {sqrt(\$sum2/\$n - (\$sum/\$n)**2)}]
}
```

}}

1. Demonstration

foreach val {2 4 4 4 5 5 7 9} {

```   set sd [sd \$val]
```

} sd stop puts "the standard deviation is: \$sd"</lang>

XPL0

<lang XPL0>include c:\cxpl\codes; \intrinsic 'code' declarations int A, I; real N, S, S2; [A:= [2,4,4,4,5,5,7,9]; N:= 0.0; S:= 0.0; S2:= 0.0; for I:= 0 to 8-1 do

```       [N:= N + 1.0;
S:= S + float(A(I));
S2:= S2 + float(sq(A(I)));
RlOut(0, sqrt((S2/N) - sq(S/N)));
];
```

CrLf(0); ]</lang>

Output:

```    0.00000    1.00000    0.94281    0.86603    0.97980    1.00000    1.39971    2.00000
```