Standard deviation

From Rosetta Code
Jump to: navigation, search
Task
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, {2,4,4,4,5,5,7,9}, which is 2.

See also: Moving Average

Contents

[edit] 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;

Sample output:

Deviation 2.00000E+00

[edit] 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.

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

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

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

[edit] AutoHotkey

ahk forum: discussion

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

[edit] AWK

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

output:

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

[edit] Axiom

We implement a domain with dependent type T with the operation + and identity 0:
)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)
This can be called using:
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)

[edit] BBC BASIC

Uses the MOD(array()) and SUM(array()) functions.

      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)

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

[edit] C

#include <stdio.h>
#include <stdlib.h>
#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;
}
#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);
}
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;
}

[edit] C++

#include <algorithm>
#include <iostream>
#include <iterator>
#include <cmath>
#include <vector>
#include <iterator>
#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 ;
}

[edit] C#

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);
}
}
}
0
1
0,942809041582063
0,866025403784439
0,979795897113271
1
1,39970842444753
2

[edit] Clojure

 
(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
 
 

[edit] COBOL

Using an intrinsic function:

FUNCTION STANDARD-DEVIATION(2, 4, 4, 4, 5, 5, 7, 9)

How this is implemented in the standard:

FUNCTION SQRT(FUNCTION VARIANCE(2, 4, 4, 4, 5, 5, 7, 9))

A complete implementation:

Works with: OpenCOBOL version 2.0
       >>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.

[edit] CoffeeScript

Uses a class instance to maintain state.

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

 

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

[edit] Common 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)))
 

Based on some googled web site; written ages ago.

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

[edit] Component Pascal

BlackBox Component Builder

 
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.
 

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

[edit] 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());
}
}

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

[edit] 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:

program StdDevTest;
 
{$APPTYPE CONSOLE}
 
uses
Math;
 
begin
Writeln(PopnStdDev([2,4,4,4,5,5,7,9]));
Readln;
end.

[edit] 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.
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]
}
? def [insert, stddev] := makeRunningStdDev()
# value: <insert>, <stddev>
 
? [stddev()]
# 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

[edit] Emacs Lisp

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

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

[edit] 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 ).
 
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

[edit] 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 ;

[edit] 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 ;

This variation is more numerically stable when there are large numbers of samples or large sample ranges.

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

[edit] 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.

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

[edit] 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
 
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
 

[edit] Go

Algorithm to reduce rounding errors from WP article.

State maintained with a closure.

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

Output:

0
1
0.9428090415820634
0.8660254037844386
0.9797958971132713
1
1.3997084244475304
2

[edit] Groovy

Solution:

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

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

[edit] Haskell

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

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]

[edit] 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));
}
}
0
1
0.942809041582063
0.866025403784439
0.979795897113271
1
1.39970842444753
2

[edit] 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
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

[edit] Icon and Unicon

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

[edit] 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.

   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

Alternatives:
Using verbose names for J primitives.

   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
Translation of: R

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

   require'stats'
(%:@:(%~<:)@:# * stddev)\ 2 4 4 4 5 5 7 9
0 1 0.942809 0.866025 0.979796 1 1.39971 2

[edit] 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));
}
}
}

[edit] JavaScript

Uses a closure.

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(', ');

output:

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

[edit] jq

[edit] 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.
# Compute the standard deviation of the observations 
# 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 };
 
# 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);
 
# Begin:
simulate

Example 1

# observations.txt
2
4
4
4
5
5
7
9
Output:
 
$ jq -s -f Dynamic_standard_deviation.jq observations.txt
0
1
0.9428090415820634
0.8660254037844386
0.9797958971132711
0.9999999999999999
1.3997084244475302
1.9999999999999998
 

[edit] 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.

# requires jq version > 1.4
def simulate(stream):
foreach stream as $observation
(initial_state;
update_state($observation);
standard_deviation);

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

[edit] 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.

#!/bin/bash
 
# 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 };
 
# 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

Example 3

$ ./standard_deviation_server.sh
Next observation: 10
0
Next observation: 20
5
Next observation: 0
8.16496580927726
 

[edit] Liberty BASIC

Using a global array to maintain the state. Implements definition explicitly.

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

[edit] Lua

Uses a closure. Translation of JavaScript.

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

[edit] Mathematica

runningSTDDev[n_] := (If[Not[ValueQ[$Data]], $Data = {}];
StandardDeviation[AppendTo[$Data, n]])


[edit] MATLAB / Octave

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

  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

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

  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
 
 

[edit] МК-61/52

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

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

[edit] Nimrod

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)

Output:

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

[edit] Objective-C

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

[edit] Blocks

Works with: Mac OS X version 10.6+
Works with: iOS version 4+
#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;
}

[edit] Objeck

Translation of: Java
 
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();
}
}
}
 

[edit] 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)

Sample output:

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

[edit] ooRexx

Works with: oorexx
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

[edit] PARI/GP

Uses the Cramer-Young updating algorithm. For demonstration it displays the mean and variance at each step.

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

[edit] 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;
}
}
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";

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

# <(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);
Output:
0
1
0.942809041582063
0.866025403784439
0.979795897113272
1
1.39970842444753
2

[edit] Perl 6

Works with: Rakudo Star version 2010.08

Using a closure:

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;

Using a state variable:

# 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>;
Output:
0
1
0.942809041582063
0.866025403784439
0.979795897113271
1
1.39970842444753
2

[edit] PL/I

*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;
Output:
AVERAGE= 5.50000E+0000;
Standard deviation       2.87228E+0000 

[edit] 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)) ) )

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

[edit] PowerShell

This implementation takes the form of an advanced function which can act like a cmdlet and receive input from the pipeline.

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

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

[edit] 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

Outputs

0.0000000000
1.0000000000
0.9428090416
0.8660254038
0.9797958971
1.0000000000
1.3997084244
2.0000000000

[edit] Python

[edit] 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

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

[edit] Using a class instance

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

[edit] Using a Closure

Works with: Python version 3.x
>>> 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

[edit] Using an extended generator

Works with: Python version 2.5+
>>> 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

[edit] R

[edit] Built-in Std Dev fn

#The built-in standard deviation function applies the Bessel correction.  To reverse this, we can apply an uncorrection.
#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

[edit] From scratch

#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

[edit] Racket

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

[edit] REXX

Uses running sums.

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

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

[edit] Ruby

[edit] 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.

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

[edit] Closure

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), ", "}
0.0, 1.0, 0.942809041582063, 0.866025403784439, 0.979795897113272, 1.0, 1.39970842444753, 2.0, 


[edit] Run BASIC

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


[edit] Scala

Library: 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))))
}

[edit] 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;
 
 

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

[edit] 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))))
 

[edit] Smalltalk

Works with: GNU 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 ]
].
|sdacc sd|
sdacc := SDAccum new.
 
#( 2 4 4 4 5 5 7 9 ) do: [ :v | sd := sdacc value: v ].
('std dev = %1' % { sd }) displayNl.

[edit] 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()

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

[edit] Tcl

[edit] With a Class

Works with: Tcl version 8.6
or
Library: TclOO
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])}
}
}
 
# 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"

which produces the output:

the standard deviation is: 2.0

[edit] With a Coroutine

Works with: Tcl version 8.6
# 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)}]
}
}}
 
# Demonstration
foreach val {2 4 4 4 5 5 7 9} {
set sd [sd $val]
}
sd stop
puts "the standard deviation is: $sd"

[edit] Visual Basic

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

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

Output:

0
1
0.942809041582063
0.866025403784439
0.979795897113271
1
1.39970842444753
2

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

Output:

    0.00000    1.00000    0.94281    0.86603    0.97980    1.00000    1.39971    2.00000

[edit] 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())
}
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
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox