Cumulative standard deviation: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎Using a class instance: Add callable class)
m (→‎{{header|Python}}: Distinguish sub-headers.)
Line 2,381: Line 2,381:


=={{header|Python}}==
=={{header|Python}}==
===Using a function with attached properties===
===Python: Using a function with attached properties===
The program should work with Python 2.x and 3.x,
The program should work with Python 2.x and 3.x,
although the output would not be a tuple in 3.x
although the output would not be a tuple in 3.x
Line 2,425: Line 2,425:
You could rename the method <code>sd</code> to <code>__call__</code> this would make the class instance callable like a function so instead of using <code>sd_inst.sd(value)</code> it would change to <code>sd_inst(value)</code> for the same results.
You could rename the method <code>sd</code> to <code>__call__</code> this would make the class instance callable like a function so instead of using <code>sd_inst.sd(value)</code> it would change to <code>sd_inst(value)</code> for the same results.


===Using a Closure===
===Python: Using a Closure===
{{Works with|Python|3.x}}
{{Works with|Python|3.x}}
<lang python>>>> from math import sqrt
<lang python>>>> from math import sqrt
Line 2,453: Line 2,453:
9 2.0</lang>
9 2.0</lang>


===Using an extended generator===
===Python: Using an extended generator===
{{Works with|Python|2.5+}}
{{Works with|Python|2.5+}}
<lang python>>>> from math import sqrt
<lang python>>>> from math import sqrt

Revision as of 06:22, 16 June 2015

Task
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, , which is .

See also:

360 Assembly

For maximum compatibility, this program uses only the basic instruction set. Part of the code length is due to the square root algorithm and to the nice output. <lang 360asm>******** Standard deviation of a population STDDEV CSECT 

        USING  STDDEV,R13

SAVEAREA B STM-SAVEAREA(R15) 

        DC     17F'0'
        DC     CL8'STDDEV'

STM STM R14,R12,12(R13) 

        ST     R13,4(R15)
        ST     R15,8(R13)
        LR     R13,R15
        SR     R8,R8           s=0
        SR     R9,R9           ss=0
        SR     R4,R4           i=0
        LA     R6,1
        LH     R7,N

LOOPI BXH R4,R6,ENDLOOPI 

        LR     R1,R4           i
        BCTR   R1,0
        SLA    R1,1
        LH     R5,T(R1)
        ST     R5,WW           ww=t(i)
        MH     R5,=H'1000'     w=ww*1000
        AR     R8,R5           s=s+w
        LR     R15,R5
        MR     R14,R5          w*w
        AR     R9,R15          ss=ss+w*w
        LR     R14,R8          s
        SRDA   R14,32
        DR     R14,R4          /i
        ST     R15,AVG         avg=s/i
        LR     R14,R9          ss
        SRDA   R14,32
        DR     R14,R4          ss/i
        LR     R2,R15          ss/i
        LR     R15,R8          s
        MR     R14,R8          s*s
        LR     R3,R15
        LR     R15,R4          i
        MR     R14,R4          i*i
        LR     R1,R15
        LA     R14,0
        LR     R15,R3
        DR     R14,R1          (s*s)/(i*i)
        SR     R2,R15
        LR     R10,R2          std=ss/i-(s*s)/(i*i)
        LR     R11,R10         std
        SRA    R11,1           x=std/2
        LR     R12,R10         px=std

LOOPWHIL EQU * 

        CR     R12,R11         while px<>=x
        BE     ENDWHILE
        LR     R12,R11         px=x
        LR     R15,R10         std
        LA     R14,0
        DR     R14,R12         /px
        LR     R1,R12          px
        AR     R1,R15          px+std/px
        SRA    R1,1            /2
        LR     R11,R1          x=(px+std/px)/2
        B      LOOPWHIL

ENDWHILE EQU * 

        LR     R10,R11
        CVD    R4,P8           i
        MVC    C17,MASK17
        ED     C17,P8
        MVC    BUF+2(1),C17+15
        L      R1,WW
        CVD    R1,P8
        MVC    C17,MASK17
        ED     C17,P8
        MVC    BUF+10(1),C17+15
        L      R1,AVG
        CVD    R1,P8
        MVC    C18,MASK18
        ED     C18,P8
        MVC    BUF+17(5),C18+12
        CVD    R10,P8          std
        MVC    C18,MASK18
        ED     C18,P8
        MVC    BUF+31(5),C18+12
        WTO    MF=(E,WTOMSG)		  
        B      LOOPI

ENDLOOPI EQU * 

        L      R13,4(0,R13)
        LM     R14,R12,12(R13)
        XR     R15,R15
        BR     R14
        DS     0D

N DC H'8' T DC H'2',H'4',H'4',H'4',H'5',H'5',H'7',H'9' WW DS F AVG DS F P8 DS PL8 MASK17 DC C' ',13X'20',X'2120',C'-' MASK18 DC C' ',10X'20',X'2120',C'.',3X'20',C'-' C17 DS CL17 C18 DS CL18 WTOMSG DS 0F 

        DC     H'80',XL2'0000'

BUF DC CL80'N=1 ITEM=1 AVG=1.234 STDDEV=1.234 ' 

        YREGS  
        END    STDDEV</lang>
Output:
N=1  ITEM=2  AVG=2.000  STDDEV=0.000
N=2  ITEM=4  AVG=3.000  STDDEV=1.000
N=3  ITEM=4  AVG=3.333  STDDEV=0.942
N=4  ITEM=4  AVG=3.500  STDDEV=0.866
N=5  ITEM=5  AVG=3.800  STDDEV=0.979
N=6  ITEM=5  AVG=4.000  STDDEV=1.000
N=7  ITEM=7  AVG=4.428  STDDEV=1.399
N=8  ITEM=9  AVG=5.000  STDDEV=2.000

Ada

<lang ada>with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; with Ada.Text_IO; use Ada.Text_IO;

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;
  end Add;
  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
     Add (Data, Test (Item));
  end loop;
  Put_Line ("Deviation" & Float'Image (Deviation (Data)));

end Test_Deviation;</lang>

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>

AWK

<lang AWK>

  1. syntax: GAWK -f STANDARD_DEVIATION.AWK

BEGIN { 

   n = split("2,4,4,4,5,5,7,9",arr,",")
   for (i=1; i<=n; i++) {
     temp[i] = arr[i]
     printf("%g %g\n",arr[i],stdev(temp))
   }
   exit(0)

} function stdev(arr, i,n,s1,s2,variance,x) { 

   for (i in arr) {
     n++
     x = arr[i]
     s1 += x ^ 2
     s2 += x
   }
   variance = ((n * s1) - (s2 ^ 2)) / (n ^ 2)
   return(sqrt(variance))

} </lang>

Output:
2 0
4 1
4 0.942809
4 0.866025
5 0.979796
5 1
7 1.39971
9 2

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
 Implementation == R add
   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(last %,x=0) 

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

BBC BASIC

Works with: BBC BASIC for Windows

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

     FOR i% = 1 TO 8
       READ n
       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>

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

typedef struct stat_obj_struct { 

  double sum, sum2;
  size_t num;
  Action action;

} sStatObject, *StatObject;

StatObject NewStatObject( Action action ) { 

 StatObject so;

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

COBOL

Using an intrinsic function: <lang cobol>FUNCTION STANDARD-DEVIATION(2, 4, 4, 4, 5, 5, 7, 9)</lang>

How this is implemented in the standard: <lang cobol>FUNCTION SQRT(FUNCTION VARIANCE(2, 4, 4, 4, 5, 5, 7, 9))</lang>

A complete implementation:

Works with: OpenCOBOL version 2.0

<lang cobol> >>SOURCE FREE IDENTIFICATION DIVISION. PROGRAM-ID. std-dev.

ENVIRONMENT DIVISION. CONFIGURATION SECTION. REPOSITORY. 

   FUNCTION sum-arr
   .

DATA DIVISION. WORKING-STORAGE SECTION. 78 Arr-Len VALUE 8. 01 arr-area VALUE "0204040405050709". 

   03  arr                   PIC 99 OCCURS Arr-Len TIMES.

01 i PIC 99.

01 avg PIC 9(3)V99.

01 std-dev PIC 9(3)V99.

PROCEDURE DIVISION. 

   DIVIDE FUNCTION sum-arr(arr-area) BY Arr-Len GIVING avg ROUNDED

   PERFORM VARYING i FROM 1 BY 1 UNTIL i > Arr-Len
       COMPUTE arr (i) = (arr (i) - avg) ** 2
   END-PERFORM

   COMPUTE std-dev = FUNCTION SQRT(FUNCTION sum-arr(arr-area) / Arr-Len)
   DISPLAY std-dev
   .

END PROGRAM std-dev.


IDENTIFICATION DIVISION. FUNCTION-ID. sum-arr.

DATA DIVISION. LOCAL-STORAGE SECTION. 01 i PIC 99.

LINKAGE SECTION. 78 Arr-Len VALUE 8. 01 arr-area. 

   03  arr                   PIC 99 OCCURS Arr-Len TIMES.

01 arr-sum PIC 99.

PROCEDURE DIVISION USING arr-area RETURNING arr-sum. 

   INITIALIZE arr-sum *> Without this, arr-sum is initialised incorrectly.
   PERFORM VARYING i FROM 1 BY 1 UNTIL i > Arr-Len
       ADD arr (i) TO arr-sum
   END-PERFORM
   .

END FUNCTION sum-arr.</lang>

CoffeeScript

Uses a class instance to maintain state.

<lang coffeescript> class StandardDeviation 

   constructor: ->
       @sum = 0
       @sumOfSquares = 0
       @values = 0
       @deviation = 0

   include: ( n ) ->
       @values += 1
       @sum += n
       @sumOfSquares += n * n
       mean = @sum / @values
       mean *= mean
       @deviation = Math.sqrt @sumOfSquares / @values - mean

dev = new StandardDeviation values = [ 2, 4, 4, 4, 5, 5, 7, 9 ] tmp = []

for value in values 

   tmp.push value
   dev.include value
   console.log """
       Values: #{ tmp }
       Standard deviation: #{ dev.deviation }

   """

</lang>

Output:
Values: 2
Standard deviation: 0

Values: 2,4
Standard deviation: 1

Values: 2,4,4
Standard deviation: 0.9428090415820626

Values: 2,4,4,4
Standard deviation: 0.8660254037844386

Values: 2,4,4,4,5
Standard deviation: 0.9797958971132716

Values: 2,4,4,4,5,5
Standard deviation: 1

Values: 2,4,4,4,5,5,7
Standard deviation: 1.3997084244475297

Values: 2,4,4,4,5,5,7,9
Standard deviation: 2

Common Lisp

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

 (let* ((n (length samples))

(mean (/ (reduce #'+ samples) n)) (tmp (mapcar (lambda (x) (expt (- x mean) 2)) 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>

Component Pascal

BlackBox Component Builder <lang oberon2> MODULE StandardDeviation; IMPORT StdLog, Args,Strings,Math;

PROCEDURE Mean(x: ARRAY OF REAL; n: INTEGER; OUT mean: REAL); VAR i: INTEGER; total: REAL; BEGIN total := 0.0; FOR i := 0 TO n - 1 DO total := total + x[i] END; mean := total /n END Mean;

PROCEDURE SDeviation(x : ARRAY OF REAL;n: INTEGER): REAL; VAR i: INTEGER; mean,sum: REAL; BEGIN Mean(x,n,mean); sum := 0.0; FOR i := 0 TO n - 1 DO sum:= sum + ((x[i] - mean) * (x[i] - mean)); END; RETURN Math.Sqrt(sum/n); END SDeviation;

PROCEDURE Do*; VAR p: Args.Params; x: POINTER TO ARRAY OF REAL; i,done: INTEGER; BEGIN Args.Get(p); IF p.argc > 0 THEN NEW(x,p.argc); FOR i := 0 TO p.argc - 1 DO x[i] := 0.0 END; FOR i  := 0 TO p.argc - 1 DO Strings.StringToReal(p.args[i],x[i],done); StdLog.Int(i + 1);StdLog.String(" :> ");StdLog.Real(SDeviation(x,i + 1));StdLog.Ln END END END Do; END StandardDeviation. </lang> Execute: ^Q StandardDeviation.Do 2 4 4 4 5 5 7 9 ~

Output:
 1 :>  0.0
 2 :>  1.0
 3 :>  0.9428090415820634
 4 :>  0.8660254037844386
 5 :>  0.9797958971132712
 6 :>  1.0
 7 :>  1.39970842444753
 8 :>  2.0

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]) {
       stdev.addNumber(el);
       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]));
 Readln;

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>

Emacs Lisp

This implementation uses a temporary buffer (the central data structure of emacs) to have simple local variables.

<lang lisp>(defun running-std (x) 

 ; ensure that we have a float to avoid potential integer math errors.
 (setq x (float x))
 ; define variables to use
 (defvar running-sum 0 "the running sum of all known values")
 (defvar running-len 0 "the running number of all known values")
 (defvar running-squared-sum 0 "the running squared sum of all known values")
 ; and make them local to this buffer
 (make-local-variable 'running-sum)
 (make-local-variable 'running-len)
 (make-local-variable 'running-squared-sum)
 ; now process the new value
 (setq running-sum (+ running-sum x))
 (setq running-len (1+ running-len))
 (setq running-squared-sum (+ running-squared-sum (* x x)))
 ; and calculate the new standard deviation
 (sqrt (- (/ running-squared-sum 
             running-len) (/ (* running-sum running-sum) 
                                (* running-len running-len )))))</lang>

<lang lisp>(with-temp-buffer 

 (loop for i in '(2 4 4 4 5 5 7 9) do 
       (insert (number-to-string (running-std i)))
       (newline))
 (message (buffer-substring (point-min) (1- (point-max)))))

"0.0 1.0 0.9428090415820636 0.8660254037844386 0.9797958971132716 1.0 1.399708424447531 2.0"</lang>

Erlang

<lang Erlang> -module( standard_deviation ).

-export( [add_sample/2, create/0, destroy/1, get/1, task/0] ).

-compile({no_auto_import,[get/1]}).

add_sample( Pid, N ) -> Pid ! {add, N}.

create() -> erlang:spawn_link( fun() -> loop( [] ) end ).

destroy( Pid ) -> Pid ! stop.

get( Pid ) -> Pid ! {get, erlang:self()}, receive {get, Value, Pid} -> Value end.

task() -> Pid = create(), [add_print(Pid, X, add_sample(Pid, X)) || X <- [2,4,4,4,5,5,7,9]], destroy( Pid ).


add_print( Pid, N, _Add ) -> io:fwrite( "Standard deviation ~p when adding ~p~n", [get(Pid), N] ).

loop( Ns ) -> receive {add, N} -> loop( [N | Ns] ); {get, Pid} -> Pid ! {get, loop_calculate( Ns ), erlang:self()}, loop( Ns ); stop -> ok end.

loop_calculate( Ns ) -> Average = loop_calculate_average( Ns ), math:sqrt( loop_calculate_average([math:pow(X - Average, 2) || X <- Ns]) ).

loop_calculate_average( Ns ) -> lists:sum( Ns ) / erlang:length( Ns ). </lang>

Output:
9> standard_deviation:task().
Standard deviation 0.0 when adding 2
Standard deviation 1.0 when adding 4
Standard deviation 0.9428090415820634 when adding 4
Standard deviation 0.8660254037844386 when adding 4
Standard deviation 0.9797958971132712 when adding 5
Standard deviation 1.0 when adding 5
Standard deviation 1.3997084244475302 when adding 7
Standard deviation 2.0 when adding 9

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,

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

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>

Old style, two ways

Early computers loaded the entire programme and its working storage into memory and left it there throughout the run. Uninitialised variables would start with whatever had been left in memory at their address by whatever last used those addresses, though some systems would clear all of memory to zero or possibly some other value before each load. Either way, if a routine was invoked a second time, its variables would have the values left in them by their previous invocation. The DATA statement allows initial values to be specified, and allows repeat counts as well. It is not an executable statement: it is not re-executed on second and subsequent invocations of the containing routine. Thus, it is easy to have a routine employ counters and the like, visible only within themselves and initialised to zero or whatever suited.

With more complex operating systems, routines that relied on retaining values across invocations might no longer work - perhaps a fresh version of the routine would be loaded to memory (perhaps at odd intervals), or, on exit, the working storage would be discarded. There was a half-way scheme, whereby variables that had appeared in DATA statements would be retained while the others would be discarded. This subtle indication has been discarded in favour of the explicit SAVE statement, naming those variables whose values are to be retained between invocations, though compilers might also offer an option such as "automatic" (for each invocation always allocate then discard working memory) and "static" (retain values), possibly introducing non-standard keywords as well. Otherwise, the routines would have to use storage global to them such as additional parameters, or, COMMON storage and in later Fortran, the MODULE arrangements for shared items. The persistence of such storage can still be limited, but by naming them in the main line can be ensured for the life of the run. The other routines with access to such storage could enable re-initialisation, additional reports, or multiple accumulations, etc.

Since the standard deviation can be calculated in a single pass through the data, producing values for the standard deviation of all values so far supplied is easily done without re-calculation. Accuracy is quite another matter. Calculations using deviances from a working mean are much better, and capturing the first X as the working mean would be easy, just test on N = 0. The sum and sum-of-squares method is quite capable of generating a negative variance, but the second method cannot, because the terms going added in to V are never negative. This could be explored by comparing the results computed from StdDev(A), StdDev(A + 10), StdDev(A + 100), StdDev(A + 1000), etc.

Incidentally, Fortran implementations rarely enable reentrancy for the WRITE statement, so, since here the functions are invoked in a WRITE statement, the functions cannot themselves use WRITE statements, say for debugging. <lang Fortran> 

     REAL FUNCTION STDDEV(X)	!Standard deviation for successive values.
      REAL X		!The latest value.
      REAL V		!Scratchpad.
      INTEGER N	!Ongoing: count of the values.
      REAL EX,EX2	!Ongoing: sum of X and X**2.
      SAVE N,EX,EX2		!Retain values from one invocation to the next.
      DATA N,EX,EX2/0,0.0,0.0/	!Initial values.
       N = N + 1		!Another value arrives.
       EX = X + EX		!Augment the total.
       EX2 = X**2 + EX2	!Augment the sum of squares.
       STDDEV = SQRT(EX2/N - (EX/N)**2)	!The variance, but, it might come out negative!
     END FUNCTION STDDEV	!For the sequence of received X values.

     REAL FUNCTION STDDEVP(X)	!Standard deviation for successive values.
      REAL X		!The latest value.
      INTEGER N	!Ongoing: count of the values.
      REAL A,V		!Ongoing: average, and sum of squared deviations.
      SAVE N,A,V		!Retain values from one invocation to the next.
      DATA N,A,V/0,0.0,0.0/	!Initial values.
       N = N + 1			!Another value arrives.
       V = (N - 1)*(X - A)**2 /N + V	!First, as it requires the existing average.
       A = (X - A)/N + A		!= [x + (n - 1).A)]/n: recover the total from the average.
       STDDEVP = SQRT(V/N)	!V can never be negative, even with limited precision.
     END FUNCTION STDDEVP	!For the sequence of received X values.

     PROGRAM TEST
     INTEGER I		!A stepper.
     REAL A(8)		!The example data.
     DATA A/2.0,3*4.0,2*5.0,7.0,9.0/	!Alas, another opportunity to use @ passed over.
     WRITE (6,1)
   1 FORMAT ("Progressive calculation of the standard deviation."/
    1 " I, A(I),   EX EX2,   Av V*N.")
     DO I = 1,8			!Step along the data series,
       WRITE (6,2) I,A(I),STDDEV(A(I)),STDDEVP(A(I))	!Showing progressive values.
   2   FORMAT (I2,F6.1,2F10.6)		!Should do for the example.
     END DO				!On to the next value.
     END

</lang>

Output: 

Progressive calculation of the standard deviation.
 I, A(I),   EX EX2,   Av V*N.
 1   2.0  0.000000  0.000000
 2   4.0  1.000000  1.000000
 3   4.0  0.942809  0.942809
 4   4.0  0.866025  0.866025
 5   5.0  0.979796  0.979796
 6   5.0  1.000000  1.000000
 7   7.0  1.399708  1.399708
 8   9.0  2.000000  2.000000

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

Haskell

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

<lang haskell>import Data.List (genericLength) import Data.STRef import Control.Monad.ST

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

 where n = genericLength l
       mean = sum l / n

sdAccum :: RealFloat a => ST s (a -> ST s a) sdAccum = do 

   accum <- newSTRef []
   return $ \x -> do
       modifySTRef accum (x:)
       list <- readSTRef accum
       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
Adding 4 stdev = 1
Adding 4 stdev = 0.9428090416
Adding 4 stdev = 0.8660254038
Adding 5 stdev = 0.9797958971
Adding 5 stdev = 1
Adding 7 stdev = 1.399708424
Adding 9 stdev = 2

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>

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

<lang java>public class StdDev { 

   int n = 0;
   double sum = 0;
   double sum2 = 0;

   public double sd(double x) {

n++; sum += x; sum2 += x*x;

return Math.sqrt(sum2/n - sum*sum/n/n); 

   }

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

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

}</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

jq

Observations from a file or array

We first define a filter, "simulate", that, if given a file of observations, will emit the standard deviations of the observations seen so far. The current state is stored in a JSON object, with this structure: 

{ "n": _, "ssd": _, "mean": _ }

where "n" is the number of observations seen, "mean" is their average, and "ssd" is the sum of squared deviations about that mean.

The challenge here is to ensure accuracy for very large n, without sacrificing efficiency. The key idea in that regard is that if m is the mean of a series of n observations, x, we then have for any a: 

SIGMA( (x - a)^2 ) == SIGMA( (x-m)^2 ) + n * (a-m)^2 == SSD + n*(a-m)^2
where SSD is the sum of squared deviations about the mean.

<lang jq># Compute the standard deviation of the observations

  1. seen so far, given the current state as input:

def standard_deviation: .ssd / .n | sqrt;

def update_state(observation): 

 def sq: .*.;
 ((.mean * .n + observation) / (.n + 1)) as $newmean
 | (.ssd + .n * ((.mean - $newmean) | sq)) as $ssd
 | { "n": (.n + 1),
     "ssd":  ($ssd + ((observation - $newmean) | sq)),
     "mean": $newmean }

def initial_state: { "n": 0, "ssd": 0, "mean": 0 };

  1. Given an array of observations presented as input:

def simulate: 

 def _simulate(i; observations):
   if (observations|length) <= i then empty
   else update_state(observations[i]) 
     | standard_deviation, _simulate(i+1; observations)
   end ;
 . as $in | initial_state | _simulate(0; $in);
  1. Begin:

simulate</lang> Example 1 

# observations.txt
2
4
4
4
5
5
7
9
Output:

<lang sh> $ jq -s -f Dynamic_standard_deviation.jq observations.txt 0 1 0.9428090415820634 0.8660254037844386 0.9797958971132711 0.9999999999999999 1.3997084244475302 1.9999999999999998 </lang>

Observations from a stream

Recent versions of jq (after 1.4) support retention of state while processing a stream. This means that any generator (including generators that produce items indefinitely) can be used as the source of observations, without first having to capture all the observations, e.g. in a file or array. <lang jq># requires jq version > 1.4 def simulate(stream): 

 foreach stream as $observation
   (initial_state;
    update_state($observation);
    standard_deviation);</lang>

Example 2: 

simulate( range(0;10) )
Output:
0
0.5
0.816496580927726
1.118033988749895
1.4142135623730951
1.707825127659933
2
2.29128784747792
2.581988897471611
2.8722813232690143

Observations from an external stream

The following illustrates how jq can be used to process observations from an external (potentially unbounded) stream, one at a time. Here we use bash to manage the calls to jq.

The definitions of the filters update_state/1 and initial_state/0 are as above but are repeated so that this script is self-contained. <lang sh>#!/bin/bash

  1. jq is assumed to be on PATH

PROGRAM=' def standard_deviation: .ssd / .n | sqrt;

def update_state(observation): 

 def sq: .*.;
 ((.mean * .n + observation) / (.n + 1)) as $newmean
 | (.ssd + .n * ((.mean - $newmean) | sq)) as $ssd
 | { "n": (.n + 1),
     "ssd":  ($ssd + ((observation - $newmean) | sq)),
     "mean": $newmean }

def initial_state: { "n": 0, "ssd": 0, "mean": 0 };

  1. Input should be [observation, null] or [observation, state]

def standard_deviations: 

 . as $in
 | if type == "array" then
     (if .[1] == null then initial_state else .[1] end) as $state
     | $state | update_state($in[0])
     | standard_deviation, .
   else empty
   end

standard_deviations ' state=null while read -p "Next observation: " observation do 

 result=$(echo "[ $observation, $state ]" | jq -c "$PROGRAM")
 sed -n 1p <<< "$result"
 state=$(sed -n 2p <<< "$result")

done</lang> Example 3 <lang sh>$ ./standard_deviation_server.sh Next observation: 10 0 Next observation: 20 5 Next observation: 0 8.16496580927726 </lang>

Julia

Use a closure to create a running standard deviation function. <lang Julia> function makerunningstd() 

   a = zero(Float64)
   b = zero(Float64)
   n = zero(Int64)
   function runningstd(x)
       a += x
       b += x^2
       n += 1
       std = sqrt(n*b - a^2)/n
       return std
   end
   return runningstd

end

test = [2, 4, 4, 4, 5, 5, 7, 9]

rstd = makerunningstd()

println("Perform a running standard deviation of ", test) for i in test 

   println(i, " => ", rstd(i))

end </lang>

Output:
Perform a running standard deviation of [2,4,4,4,5,5,7,9]
2 => 0.0
4 => 1.0
4 => 0.9428090415820635
4 => 0.8660254037844386
5 => 0.9797958971132712
5 => 1.0
7 => 1.3997084244475302
9 => 2.0

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
       read x
       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>

Here is a vectorized one line solution as a function <lang Matlab> function stdDevEval(n) disp(sqrt(sum((n-sum(n)/length(n)).^2)/length(n))); end </lang>

МК-61/52

<lang>0 П4 П5 П6 С/П П0 ИП5 + П5 ИП0 x^2 ИП6 + П6 КИП4 ИП6 ИП4 / ИП5 ИП4 / x^2 - КвКор БП 04</lang>

Instruction: В/О С/П number С/П number С/П ...

Nim

<lang nim>import math, strutils

var sdSum, sdSum2, sdN = 0.0 proc sd(x): float = 

 sdN    += 1
 sdSum  += float(x)
 sdSum2 += float(x*x)
 sqrt(sdSum2/sdN - sdSum*sdSum/sdN/sdN)

for value in [2,4,4,4,5,5,7,9]: 

 echo value, " ", formatFloat(sd(value), precision = 0)</lang>
Output:
2 0
4 1
4 0.942809
4 0.866025
5 0.979796
5 1
7 1.39971
9 2

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() { 

 @autoreleasepool {

   double v[] = { 2,4,4,4,5,5,7,9 };

   SDAccum *sdacc = [[SDAccum alloc] init];

   for(int i=0; i < sizeof(v)/sizeof(*v) ; i++)
     printf("adding %f\tstddev = %f\n", v[i], [sdacc value: v[i]]);

 }
 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);
 };

}

int main() { 

 @autoreleasepool {

   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]));

 }
 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->AddNum(test_data[i]);
       sd->GetSD()->PrintLine();
     };
   }
   
   method : public : AddNum(num : Float) ~ Nil {
     nums->AddBack(num);
   }
   
   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>
Output:
List: 2  4  4  4  5  5  7  9
Standard deviation: 2

ooRexx

Works with: oorexx

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

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

}; addpoint(x)={ 

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

}; addpoints(v)={ 

 print(newpoint(v[1]));
 for(i=2,#v,print(addpoint(v[i])));
 print("Mean: ",myT/myN);
 print("Standard deviation: ",sqrt(myS/myN))

}; addpoints([2,4,4,4,5,5,7,9])</lang>

Pascal

Translation of: AWK

<lang pascal>program stddev; uses math; const 

 n=8;

var 

 arr: array[1..n] of real =(2,4,4,4,5,5,7,9);

function stddev(n: integer): real; var 

  i: integer;
  s1,s2,variance,x: real;

begin 

   for i:=1 to n do
   begin
     x:=arr[i];
     s1:=s1+power(x,2);
     s2:=s2+x
   end;
   variance:=((n*s1)-(power(s2,2)))/(power(n,2));
   stddev:=sqrt(variance)

end; var 

  i: integer;

begin 

   for i:=1 to n do
   begin
     writeln(i,' item=',arr[i]:2:0,' stddev=',stddev(i):18:15)
   end

end.</lang>

Output:
1 item= 2 stddev= 0.000000000000000
2 item= 4 stddev= 1.000000000000000
3 item= 4 stddev= 0.942809041582064
4 item= 4 stddev= 0.866025403784439
5 item= 5 stddev= 0.979795897113271
6 item= 5 stddev= 1.000000000000000
7 item= 7 stddev= 1.399708424447530
8 item= 9 stddev= 2.000000000000000

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>² { 

   my $num, $sum, $sum2;
   sub stddev {

my $x = shift; $num++; return sqrt( ($sum2 += $x**2) / $num - (($sum += $x) / $num)**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 R/ [+] @a;
   sqrt @a R/ [+] map (* - $mean)**2, @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) { 

   sqrt
       (.[2] += $x**2) / ++.[0] -
       ((.[1] += $x) / .[0])**2
   given state @;

}

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

PHP

This is just straight PHP class usage, respecting the specifications "stateful" and "one at a time": <lang PHP><?php class sdcalc { 

   private  $cnt, $sumup, $square;

   function __construct() {
      $this->reset();
   }
   # callable on an instance
   function reset() {
      $this->cnt=0; $this->sumup=0; $this->square=0;
   }
   function add($f) {
       $this->cnt++;
       $this->sumup  += $f;
       $this->square += pow($f, 2);
       return $this->calc();
   }
   function calc() {
       if ($this->cnt==0 || $this->sumup==0) {
           return 0;
       } else {
           return sqrt($this->square / $this->cnt - pow(($this->sumup / $this->cnt),2));
       }
   }
}
  1. start test, adding test data one by one

$c = new sdcalc(); foreach ([2,4,4,4,5,5,7,9] as $v) { 

   printf('Adding %g: result %g%s', $v, $c->add($v), PHP_EOL);

}</lang>

This will produce the output: 

Adding 2: result 0
Adding 4: result 1
Adding 4: result 0.942809
Adding 4: result 0.866025
Adding 5: result 0.979796
Adding 5: result 1
Adding 7: result 1.39971
Adding 9: result 2

PL/I

<lang pli>*process source attributes xref; 

stddev: proc options(main);
  declare a(10) float init(1,2,3,4,5,6,7,8,9,10);
  declare stdev float;
  declare i fixed binary;   
 
  stdev=std_dev(a);
  put skip list('Standard deviation', stdev);   
 
  std_dev: procedure(a) returns(float);
    declare a(*) float, n fixed binary;
    n=hbound(a,1);
    begin;
      declare b(n) float, average float;
      declare i fixed binary;
      do i=1 to n;
        b(i)=a(i);
      end;
      average=sum(a)/n;
      put skip data(average);
      return( sqrt(sum(b**2)/n - average**2) );
    end;
  end std_dev;

end;</lang>
Output:
AVERAGE= 5.50000E+0000;
Standard deviation       2.87228E+0000

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
   Read.i x
   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
 Static in_kvadrater.q
 in_summa+in
 in_kvadrater+in*in
 antal+1
 ProcedureReturn Pow((in_kvadrater/antal)-Pow(in_summa/antal,2),0.50)

EndProcedure

data section

DataSection MyList: 

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

EndDataSection</lang>

Output:
 0.0000000000
 1.0000000000
 0.9428090416
 0.8660254038
 0.9797958971
 1.0000000000
 1.3997084244
 2.0000000000

Python

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>

Python: 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>

Python: Callable class

You could rename the method sd to __call__ this would make the class instance callable like a function so instead of using sd_inst.sd(value) it would change to sd_inst(value) for the same results.

Python: 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>

Python: Using an extended generator

Works with: Python version 2.5+

<lang python>>>> from math import sqrt >>> def sdcreator(): sum = sum2 = n = 0 while True: x = yield sqrt(sum2/n - sum*sum/n/n) if n else None

sum += x sum2 += x*x n += 1.0

>>> sd = sdcreator() >>> sd.send(None) >>> for value in (2,4,4,4,5,5,7,9): print (value, sd.send(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 rsplus>#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 rsplus>#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>

Racket

<lang racket>

  1. lang racket

(require math) (define running-stddev 

 (let ([ns '()])
   (λ(n) (set! ns (cons n ns)) (stddev ns))))
run it on each number, return the last result

(last (map running-stddev '(2 4 4 4 5 5 7 9))) </lang>

REXX

Uses running sums. <lang rexx>/*REXX pgm finds & displays the standard deviation of a given set of #s.*/ parse arg # /*any optional args on the C.L. ?*/ if #= then #=2 4 4 4 5 5 7 9 /*None given? Then use default.*/ w=words(#); s=0; ss=0 /*define: #items; a couple sums.*/ 

           do j=1  for w;       _=word(#,j);      s=s+_;        ss=ss+_*_
           say  '   item'  right(j,length(w))":"  right(_,4),
                '   average='  left(s/j,12),
                '   standard deviation=' left(sqrt( ss/j - (s/j)**2 ),15)
           end   /*j*/

exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────SQRT subroutine─────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits 11 numeric form; m.=11; p=d+d%4+2; parse value format(x,2,1,,0) 'E0' with g 'E' _ . g=g*.5'E'_%2; do j=0 while p>9; m.j=p; p=p%2+1; end 

do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g); end

numeric digits d; return g/1</lang>

Output:

using the default input 

   item 1:    2    average= 2               standard deviation= 0
   item 2:    4    average= 3               standard deviation= 1
   item 3:    4    average= 3.33333333      standard deviation= 0.942809047
   item 4:    4    average= 3.5             standard deviation= 0.866025404
   item 5:    5    average= 3.8             standard deviation= 0.979795897
   item 6:    5    average= 4               standard deviation= 1
   item 7:    7    average= 4.42857143      standard deviation= 1.39970843
   item 8:    9    average= 5               standard deviation= 2

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,


Run BASIC

<lang runbasic>dim sdSave$(100) 'can call up to 100 versions 

                 'holds (space-separated) number of data , sum of values and sum of squares

sd$ = "2,4,4,4,5,5,7,9"

for num = 1 to 8 

stdData = val(word$(sd$,num,","))
 sumVal = sumVal + stdData
 sumSqs = sumSqs + stdData^2

 ' standard deviation = square root of (the average of the squares less the square of the average)
 standDev   =((sumSqs / num) - (sumVal /num) ^ 2) ^ 0.5

 sdSave$(num) = str$(num);" ";str$(sumVal);" ";str$(sumSqs)
 print num;" value in = ";stdData; " Stand Dev = "; using("###.######", standDev)

next num</lang> 

1 value in = 2 Stand Dev =   0.000000
2 value in = 4 Stand Dev =   1.000000
3 value in = 4 Stand Dev =   0.942809
4 value in = 4 Stand Dev =   0.866025
5 value in = 5 Stand Dev =   0.979796
6 value in = 5 Stand Dev =   1.000000
7 value in = 7 Stand Dev =   1.399708
8 value in = 9 Stand Dev =   2.000000

SAS

<lang SAS>

  • --Load the test data;

data test1; 

  input x @@;
  obs=_n_;

datalines; 2 4 4 4 5 5 7 9

run;

  • --Create a dataset with the cummulative data for each set of data for which the SD should be calculated;

data test2 (drop=i obs); 

  set test1;
  y=x;
  do i=1 to n;
     set test1 (rename=(obs=setid)) nobs=n point=i;
     if obs<=setid then output;
  end;

proc sort; 

  by setid;

run;

  • --Calulate the standards deviation (and mean) using PROC MEANS;

proc means data=test2 vardef=n noprint; *--use vardef=n option to calculate the population SD; 

  by setid;
  var y;
  output out=stat1 n=n mean=mean std=sd;

run;

  • --Output the calculated standard deviations;

proc print data=stat1 noobs; 

  var n sd /*mean*/;

run; </lang>

Output:
N       SD

1    0.00000
2    1.00000
3    0.94281
4    0.86603
5    0.97980
6    1.00000
7    1.39971
8    2.00000


Scala

Library: Scala

<lang Scala>import scala.math._ import Numeric.Implicits._

object StddevCalc extends App { 

 def avg[T](ts: Iterable[T])(implicit num: Fractional[T]): T = {
   num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T
 }

 def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = {
   val mean = avg(ts) // Leave val type of T
   val stdDev = // Root of mean diffs
     sqrt(num.toDouble(ts.foldLeft(num.zero)((b, a) => num.plus(b, num.times(num.minus(a, mean), num.minus(a, mean))))) / ts.size)
   (mean, stdDev)
 }

 def calcAvgAndStddev(ts: Iterable[BigDecimal]): (Double, Double) = // Overloaded for BigDecimal
   calcAvgAndStddev(ts.map(_.toDouble))

 println(calcAvgAndStddev(List(2.0E0, 4.0, 4, 4, 5, 5, 7, 9)))
 println(calcAvgAndStddev(Set(1.0, 2, 3, 4)))
 println(calcAvgAndStddev(0.1 to 1.1 by 0.05))
 println(calcAvgAndStddev(List(BigDecimal(120), BigDecimal(1200))))

}</lang>

Scheme

Works with: Racket
This example is incorrect. Please fix the code and remove this message.

Details: running-stddev procedure only in Racket? Usage example missing, and not evident.

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

Scilab

Scilab has the built-in function stdev to compute the standard deviation of a sample so it is straightforward to have the standard deviation of a sample with a correction of the bias. <lang>T=[2,4,4,4,5,5,7,9]; stdev(T)*sqrt((length(T)-1)/length(T))</lang>

Output:
-->T=[2,4,4,4,5,5,7,9];
-->stdev(T)*sqrt((length(T)-1)/length(T))
   ans  =     2.

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>

SQL

Works with: Postgresql

<lang SQL>-- the minimal table create table if not exists teststd (n double precision not null);

-- code modularity with view, we could have used a common table expression instead create view vteststd as 

 select count(n) as cnt,
 sum(n) as tsum,
 sum(power(n,2)) as tsqr

from teststd;

-- you can of course put this code into every query create or replace function std_dev() returns double precision as $$ 

select sqrt(tsqr/cnt - (tsum/cnt)^2) from vteststd;

$$ language sql;

-- test data is: 2,4,4,4,5,5,7,9 insert into teststd values (2); select std_dev() as std_deviation; insert into teststd values (4); select std_dev() as std_deviation; insert into teststd values (4); select std_dev() as std_deviation; insert into teststd values (4); select std_dev() as std_deviation; insert into teststd values (5); select std_dev() as std_deviation; insert into teststd values (5); select std_dev() as std_deviation; insert into teststd values (7); select std_dev() as std_deviation; insert into teststd values (9); select std_dev() as std_deviation; -- cleanup test data delete from teststd; </lang> With a command like psql <rosetta-std-dev.sql you will get an output like this: (duplicate lines generously deleted, locale is DE) 

CREATE TABLE
FEHLER:  Relation »vteststd« existiert bereits
CREATE FUNCTION
INSERT 0 1
 std_deviation 
---------------
             0
(1 Zeile)

INSERT 0 1
 std_deviation 
---------------
             1
 0.942809041582063
 0.866025403784439
 0.979795897113272
             1
 1.39970842444753
             2
DELETE 8

Swift

<lang Swift>import Darwin class stdDev{ 

   var n:Double = 0.0
   var sum:Double = 0.0
   var sum2:Double = 0.0
   
   init(){
       
       let testData:[Double] = [2,4,4,4,5,5,7,9];
       for x in testData{
           
           var a:Double = calcSd(x)
           println("value \(Int(x)) SD = \(a)");
       }
       
   }
   
   func calcSd(x:Double)->Double{
       
       n += 1
       sum += x
       sum2 += x*x
       return sqrt( sum2 / n - sum*sum / n / n)
   }

} var aa = stdDev()</lang>

Output:
value 2 SD = 0.0
value 4 SD = 1.0
value 4 SD = 0.942809041582063
value 4 SD = 0.866025403784439
value 5 SD = 0.979795897113271
value 5 SD = 1.0
value 7 SD = 1.39970842444753
value 9 SD = 2.0

Functional:

<lang Swift> func standardDeviation(arr : [Double]) -> Double { 

   let length = Double(arr.count)
   let avg = arr.reduce(0, { $0 + $1 }) / length
   let sumOfSquaredAvgDiff = arr.map { pow($0 - avg, 2.0)}.reduce(0, {$0 + $1})
   return sqrt(sumOfSquaredAvgDiff / length)

}

let responseTimes = [ 18.0, 21.0, 41.0, 42.0, 48.0, 50.0, 55.0, 90.0 ]

standardDeviation(responseTimes) // 20.8742514835862 standardDeviation([2,4,4,4,5,5,7,9]) // 2.0 </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>

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>

TI-83 BASIC

On the TI-83 family, standard deviation of a population is a builtin function (σx): 

• Press [STAT] select [EDIT] followed by [ENTER]
• then enter for list L1 in the table : 2, 4, 4, 4, 5, 5, 7, 9
• Or enter {2,4,4,4,5,5,7,9}→L1
• Press [STAT] select [CALC] then [1-Var Stats] select list L1 followed by [ENTER]
• Then σx (=2) gives the standard deviation of the population

VBScript

<lang vb>data = Array(2,4,4,4,5,5,7,9)

For i = 0 To UBound(data) WScript.StdOut.Write "value = " & data(i) &_ " running sd = " & sd(data,i) WScript.StdOut.WriteLine Next

Function sd(arr,n) mean = 0 variance = 0 For j = 0 To n mean = mean + arr(j) Next mean = mean/(n+1) For k = 0 To n variance = variance + ((arr(k)-mean)^2) Next variance = variance/(n+1) sd = FormatNumber(Sqr(variance),6) End Function</lang>

Output:
value = 2 running sd = 0.000000
value = 4 running sd = 1.000000
value = 4 running sd = 0.942809
value = 4 running sd = 0.866025
value = 5 running sd = 0.979796
value = 5 running sd = 1.000000
value = 7 running sd = 1.399708
value = 9 running sd = 2.000000

Visual Basic

Note that the helper function avg is not named average to avoid a name conflict with WorksheetFunction.Average in MS Excel.

<lang vb>Function avg(what() As Variant) As Variant 

   'treats non-numeric strings as zero
   Dim L0 As Variant, total As Variant
   For L0 = LBound(what) To UBound(what)
       If IsNumeric(what(L0)) Then total = total + what(L0)
   Next
   avg = total / (1 + UBound(what) - LBound(what))

End Function

Function standardDeviation(fp As Variant) As Variant 

   Static list() As Variant
   Dim av As Variant, tmp As Variant, L0 As Variant

   'add to sequence if numeric
   If IsNumeric(fp) Then
       On Error GoTo makeArr   'catch undimensioned list
       ReDim Preserve list(UBound(list) + 1)
       On Error GoTo 0
       list(UBound(list)) = fp
   End If

   'get average
   av = avg(list())

   'the actual work
   For L0 = 0 To UBound(list)
       tmp = tmp + ((list(L0) - av) ^ 2)
   Next
   tmp = Sqr(tmp / (UBound(list) + 1))

   standardDeviation = tmp

   Exit Function

makeArr: 

   If 9 = Err.Number Then
       ReDim list(0)
   Else
       'something's wrong
       Err.Raise Err.Number
   End If
   Resume Next

End Function

Sub tester() 

   Dim x As Variant
   x = Array(2, 4, 4, 4, 5, 5, 7, 9)
   For L0 = 0 To UBound(x)
       Debug.Print standardDeviation(x(L0))
   Next

End Sub</lang>

Output:
 0
 1
 0.942809041582063
 0.866025403784439
 0.979795897113271
 1
 1.39970842444753
 2

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

zkl

<lang zkl>fcn sdf{ fcn(x,xs){ 

     m:=xs.append(x.toFloat()).sum(0.0)/xs.len(); 
     (xs.reduce('wrap(p,x){(x-m)*(x-m) +p},0.0)/xs.len()).sqrt() 
   }.fp1(L())

}</lang>

Output:
zkl: T(2,4,4,4,5,5,7,9).pump(Void,sdf())
2

zkl: sd:=sdf(); sd(2);sd(4);sd(4);sd(4);sd(5);sd(5);sd(7);sd(9)
2
Retrieved from "https://rosettacode.org/wiki/Cumulative_standard_deviation?oldid=116498"