Write a stateful function, class, generator or co-routine that takes a series of floating point numbers,one at a time, and returns the runningstandard 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 taskmust state which is being used.
Do not applyBessel'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.
| ||||||||||||||||||||||||||||||||||||||||||||||
T SD sum = 0.0 sum2 = 0.0 n = 0.0 F ()(x) .sum += x .sum2 += x ^ 2 .n += 1.0 R sqrt(.sum2 / .n - (.sum / .n) ^ 2)V sd_inst = SD()L(value) [2, 4, 4, 4, 5, 5, 7, 9] print(value‘ ’sd_inst(value))
2 04 14 0.9428090424 0.8660254045 0.9797958975 17 1.3997084249 2
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.
******** Standard deviation of a populationSTDDEV CSECT USING STDDEV,R13SAVEAREA 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,NLOOPI 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=stdLOOPWHIL 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 LOOPWHILENDWHILE 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 LOOPIENDLOOPI EQU * L R13,4(0,R13) LM R14,R12,12(R13) XR R15,R15 BR R14 DS 0DN DC H'8'T DC H'2',H'4',H'4',H'4',H'5',H'5',H'7',H'9'WW DS FAVG DS FP8 DS PL8MASK17 DC C' ',13X'20',X'2120',C'-' MASK18 DC C' ',10X'20',X'2120',C'.',3X'20',C'-' C17 DS CL17C18 DS CL18WTOMSG DS 0F DC H'80',XL2'0000'BUF DC CL80'N=1 ITEM=1 AVG=1.234 STDDEV=1.234 ' YREGS END STDDEV
N=1 ITEM=2 AVG=2.000 STDDEV=0.000N=2 ITEM=4 AVG=3.000 STDDEV=1.000N=3 ITEM=4 AVG=3.333 STDDEV=0.942N=4 ITEM=4 AVG=3.500 STDDEV=0.866N=5 ITEM=5 AVG=3.800 STDDEV=0.979N=6 ITEM=5 AVG=4.000 STDDEV=1.000N=7 ITEM=7 AVG=4.428 STDDEV=1.399N=8 ITEM=9 AVG=5.000 STDDEV=2.000
INCLUDE "H6:REALMATH.ACT"REAL sum,sum2INT countPROC Calc(REAL POINTER x,sd) REAL tmp1,tmp2,tmp3 RealAdd(sum,x,tmp1) ;tmp1=sum+x RealAssign(tmp1,sum) ;sum=sum+x RealMult(x,x,tmp1) ;tmp1=x*x RealAdd(sum2,tmp1,tmp2) ;tmp2=sum2+x*x RealAssign(tmp2,sum2) ;sum2=sum2+x*x count==+1 IF count=0 THEN IntToReal(0,sd) ;sd=0 ELSE IntToReal(count,tmp1) RealMult(sum,sum,tmp2) ;tmp2=sum*sum RealDiv(tmp2,tmp1,tmp3) ;tmp3=sum*sum/count RealDiv(tmp3,tmp1,tmp2) ;tmp2=sum*sum/count/count RealDiv(sum2,tmp1,tmp3) ;tmp3=sum2/count RealSub(tmp3,tmp2,tmp1) ;tmp1=sum2/count-sum*sum/count/count Sqrt(tmp1,sd) ;sd=sqrt(sum2/count-sum*sum/count/count) FIRETURNPROC Main() INT ARRAY values=[2 4 4 4 5 5 7 9] INT i REAL x,sd Put(125) PutE() ;clear screen MathInit() IntToReal(0,sum) IntToReal(0,sum2) count=0 FOR i=0 TO 7 DO IntToReal(values(i),x) Calc(x,sd) Print("x=") PrintR(x) Print(" sum=") PrintR(sum) Print(" sd=") PrintRE(sd) ODRETURNScreenshot from Atari 8-bit computer
x=2 sum=2 sd=0x=4 sum=6 sd=1x=4 sum=10 sd=.942809052x=4 sum=14 sd=.86602541x=5 sum=19 sd=.979795903x=5 sum=24 sd=1x=7 sum=31 sd=1.39970843x=9 sum=40 sd=1.99999999
withAda.Numerics.Elementary_Functions;useAda.Numerics.Elementary_Functions;withAda.Numerics.Elementary_Functions;useAda.Numerics.Elementary_Functions;withAda.Text_IO;useAda.Text_IO;withAda.Float_Text_IO;useAda.Float_Text_IO;withAda.Integer_Text_IO;useAda.Integer_Text_IO;procedureTest_DeviationistypeSampleisrecordN:Natural:=0;Sum:Float:=0.0;SumOfSquares:Float:=0.0;end record;procedureAdd(Data:inoutSample;Point:Float)isbeginData.N:=Data.N+1;Data.Sum:=Data.Sum+Point;Data.SumOfSquares:=Data.SumOfSquares+Point**2;endAdd;functionDeviation(Data:Sample)returnFloatisbeginreturnSqrt(Data.SumOfSquares/Float(Data.N)-(Data.Sum/Float(Data.N))**2);endDeviation;Data:Sample;Test:array(1..8)ofInteger:=(2,4,4,4,5,5,7,9);beginforIndexinTest'RangeloopAdd(Data,Float(Test(Index)));Put("N=");Put(Item=>Index,Width=>1);Put(" ITEM=");Put(Item=>Test(Index),Width=>1);Put(" AVG=");Put(Item=>Float(Data.Sum)/Float(Index),Fore=>1,Aft=>3,Exp=>0);Put(" STDDEV=");Put(Item=>Deviation(Data),Fore=>1,Aft=>3,Exp=>0);New_line;endloop;endTest_Deviation;
N=1 ITEM=2 AVG=2.000 STDDEV=0.000N=2 ITEM=4 AVG=3.000 STDDEV=1.000N=3 ITEM=4 AVG=3.333 STDDEV=0.943N=4 ITEM=4 AVG=3.500 STDDEV=0.866N=5 ITEM=5 AVG=3.800 STDDEV=0.980N=6 ITEM=5 AVG=4.000 STDDEV=1.000N=7 ITEM=7 AVG=4.429 STDDEV=1.400N=8 ITEM=9 AVG=5.000 STDDEV=2.000
COMMENT Cumulative standard deviation;BEGINPROCEDURE init(n, sum, sumofsquares); INTEGER n; REAL sum, sumofsquares;BEGIN n := 0; sum := 0.0; sumofsquares := 0.0END;PROCEDURE add(n, sum, sumofsquares, p); VALUE p; INTEGER n; REAL sum, sumofsquares, p;BEGIN n := n + 1; sum := sum + p; sumofsquares := sumofsquares + p * pEND;REAL PROCEDURE deviation(n, sum, sumofsquares); VALUE n, sum, sumofsquares; INTEGER n; REAL sum, sumofsquares;BEGIN REAL sumdivn; sumdivn := sum / n; deviation := sqrt(sumofsquares / n - sumdivn * sumdivn)END;REAL ARRAY testdata[1 : 8];INTEGER n, i;REAL sum, sumofsquares;testdata[1] := 2;testdata[2] := 4;testdata[3] := 4;testdata[4] := 4;testdata[5] := 5;testdata[6] := 5;testdata[7] := 7;testdata[8] := 9;init(n, sum, sumofsquares);text(1, " N ITEM AVG STDDEV");skip(1);FOR i := 1 STEP 1 UNTIL 8 DOBEGIN add(n, sum, sumofsquares, testdata[i]); rwrite(1, i, 2, 0); rwrite(1, testdata[i], 4, 0); rwrite(1, sum / i, 6, 3); rwrite(1, deviation(n, sum, sumofsquares), 6, 3); skip(1)END;END FINISH
N ITEM AVG STDDEV 1 2 2.000 0.000 2 4 3.000 1.000 3 4 3.333 0.943 4 4 3.500 0.866 5 5 3.800 0.980 6 5 4.000 1.000 7 7 4.429 1.400 8 9 5.000 2.000
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 ofALGOL 68sconformity 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); # main # ( []LONG REAL v = ( 2,4,4,4,5,5,7,9 ); LONG REAL sd; FOR i FROM LWB v TO UPB v DO sd := stat object(v[i], LOC VALUE); printf(($"value: "g(0,6)," standard dev := "g(0,6)l$, v[i], sd)) OD )
value: 2.000000 standard dev := .000000value: 4.000000 standard dev := 1.000000value: 4.000000 standard dev := .942809value: 4.000000 standard dev := .866025value: 5.000000 standard dev := .979796value: 5.000000 standard dev := 1.000000value: 7.000000 standard dev := 1.399708value: 9.000000 standard dev := 2.000000
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));# could define STDDEV as an operator for more syntatic sugar OP STDDEV = ([]LONG REAL value)LONG REAL: ( 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 ); # main # ( []LONG REAL value = ( 2,4,4,4,5,5,7,9 );# 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]; printf(($"value: "g(0,6)," standard dev := "g(0,6)l$, value[i], (stddev OF class stat)(stat))) 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)));# SKIP)
value: 2.000000 standard dev := .000000value: 4.000000 standard dev := 1.000000value: 4.000000 standard dev := .942809value: 4.000000 standard dev := .866025value: 5.000000 standard dev := .979796value: 5.000000 standard dev := 1.000000value: 7.000000 standard dev := 1.399708value: 9.000000 standard dev := 2.000000
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
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
This is an Algol W version of Oberon-07 sample which was based on the third, "unpackaged" Algol 68 sample, which was itself translated from Python.
begin record CSD ( long real sumof, sum2of, sdof; integer countof ); long real procedure sample ( reference(CSD) value sdv; long real value x) ; begin long real sum, sum2; integer n; countof(sdv) := countof(sdv) + 1; sum := sumof(sdv); sum2 := sum2of(sdv); n := countof(sdv); sum := sum + x; sum2 := sum2 + (x*x); sdof(sdv) := if n = 0 then 0 else longsqrt(sum2/n - sum*sum/n/n); sumof(sdv) := sum; sum2of(sdv) := sum2; sdof(sdv) end sample; reference(CSD) procedure NewCSD; CSD( 0, 0, 0, 0 ); begin reference(CSD) sd; sd := NewCSD; r_format := "A"; r_w := 14; r_d := 6; % set output to fixed point format % i_w := 6; s_w := 0; % also set integer and separator format % for i := 2,4,4,4,5,5,7,9 do begin long real val, newSd; val := i; newSd := sample( sd, val ); write( countof(sd), ": ", val, " -> ", newSd ) end for_i endend.
1: 2.000000 -> 0.000000 2: 4.000000 -> 1.000000 3: 4.000000 -> 0.942809 4: 4.000000 -> 0.866025 5: 5.000000 -> 0.979796 6: 5.000000 -> 1.000000 7: 7.000000 -> 1.399708 8: 9.000000 -> 2.000000
Accumulation over a fold:
-------------- CUMULATIVE STANDARD DEVIATION --------------- stdDevInc :: Accumulator -> Num -> Index -> Accumulator-- stdDevInc :: {sum:, squaresSum:, stages:} -> Real -> Integer-- -> {sum:, squaresSum:, stages:}onstdDevInc(a,n,i)setsumto(sumofa)+nsetsquaresSumto(squaresSumofa)+(n^2)setstagesto(stagesofa)&¬((squaresSum/i)-((sum/i)^2))^0.5{sum:(sumofa)+n,squaresSum:squaresSum,stages:stages}endstdDevInc--------------------------- TEST -------------------------onrunsetxsto[2,4,4,4,5,5,7,9]stagesoffoldl(stdDevInc,¬{sum:0,squaresSum:0,stages:[]},xs)--> {0.0, 1.0, 0.942809041582, 0.866025403784, 0.979795897113, 1.0, 1.399708424448, 2.0}endrun-------------------- GENERIC FUNCTIONS --------------------- foldl :: (a -> b -> a) -> a -> [b] -> aonfoldl(f,startValue,xs)tellmReturn(f)setvtostartValuesetlngtolengthofxsrepeatwithifrom1tolngsetvto|λ|(v,itemiofxs,i,xs)endrepeatreturnvendtellendfoldl-- mReturn :: First-class m => (a -> b) -> m (a -> b)onmReturn(f)-- 2nd class handler function lifted into 1st class script wrapper.ifscriptisclassoffthenfelsescriptproperty|λ|:fendscriptendifendmReturn
{0.0, 1.0, 0.942809041582, 0.866025403784, 0.979795897113, 1.0, 1.399708424448, 2.0}
Or as a map-accumulation:
-------------- CUMULATIVE STANDARD DEVIATION --------------- cumulativeStdDevns :: [Float] -> [Float]oncumulativeStdDevns(xs)scriptgoon|λ|(sq,x,i)set{s,q}tosqset_stox+sset_qtoq+(x^2){{_s,_q},((_q/i)-((_s/i)^2))^0.5}end|λ|endscriptitem2ofmapAccumL(go,{0,0},xs)endcumulativeStdDevns--------------------------- TEST -------------------------onruncumulativeStdDevns({2,4,4,4,5,5,7,9})endrun------------------------- GENERIC -------------------------- foldl :: (a -> b -> a) -> a -> [b] -> aonfoldl(f,startValue,xs)tellmReturn(f)setvtostartValuesetlngtolengthofxsrepeatwithifrom1tolngsetvto|λ|(v,itemiofxs,i,xs)endrepeatreturnvendtellendfoldl-- mapAccumL :: (acc -> x -> (acc, y)) -> acc -> [x] -> (acc, [y])onmapAccumL(f,acc,xs)-- 'The mapAccumL function behaves like a combination of map and foldl;-- it applies a function to each element of a list, passing an-- accumulating parameter from |Left| to |Right|, and returning a final-- value of this accumulator together with the new list.' (see Hoogle)scripton|λ|(a,x,i)tellmReturn(f)tosetpairto|λ|(item1ofa,x,i){item1ofpair,(item2ofa)&{item2ofpair}}end|λ|endscriptfoldl(result,{acc,[]},xs)endmapAccumL-- mReturn :: First-class m => (a -> b) -> m (a -> b)onmReturn(f)-- 2nd class handler function lifted into 1st class script wrapper.ifscriptisclassoffthenfelsescriptproperty|λ|:fendscriptendifendmReturn
{0.0, 1.0, 0.942809041582, 0.866025403784, 0.979795897113, 1.0, 1.399708424448, 2.0}arr:[]loop[24445579]'value['arr++valueprint[value"->"deviationarr]]
2 -> 0.0 4 -> 1.0 4 -> 0.9428090415820634 4 -> 0.8660254037844386 5 -> 0.9797958971132711 5 -> 0.9999999999999999 7 -> 1.39970842444753 9 -> 2.0
Data:=[2,4,4,4,5,5,7,9]fork,vinData{FileAppend,%"#"a_index" value = "v" stddev = "stddev(v)"`n",* ; send to stdout}returnstddev(x){staticn,sum,sum2n++sum+=xsum2+=x*xreturnsqrt((sum2/n)-(((sum*sum)/n)/n))}
#1 value = 2 stddev 0 0.000000#2 value = 4 stddev 0 1.000000#3 value = 4 stddev 0 0.942809#4 value = 4 stddev 0 0.866025#5 value = 5 stddev 0 0.979796#6 value = 5 stddev 0 1.000000#7 value = 7 stddev 0 1.399708#8 value = 9 stddev 0 2.000000
# syntax: GAWK -f STANDARD_DEVIATION.AWKBEGIN{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)}functionstdev(arr,i,n,s1,s2,variance,x){for(iinarr){n++x=arr[i]s1+=x^2s2+=x}variance=((n*s1)-(s2^2))/(n^2)return(sqrt(variance))}
2 04 14 0.9428094 0.8660255 0.9797965 17 1.399719 2
This example isincorrect. Please fix the code and remove this message.
|
We implement a domain with dependent type T with the operation + and identity 0:
)abbrev package TESTD TestDomainTestDomain(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 IntegerD ==> 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)
See alsoVisual Basic.
100REM Cumulative standard deviation110DECLAREEXTERNALSUBInit120DECLAREEXTERNALSUBAdd130DECLAREEXTERNALFUNCTIONDeviation140DIMTestData(1TO8)150FORI=1TO8160READTestData(I)170NEXTI180DATA2,4,4,4,5,5,7,9190CALLInit(N,Sum,SumOfSquares)200PRINT"N ITEM AVG STDDEV"210FORI=LBOUND(TestData)TOUBOUND(TestData)220CALLAdd(N,Sum,SumOfSquares,TestData(I))230PRINTUSING"# ":I;240PRINTUSING" # ":TestData(I);250PRINTUSING"#.### ":Sum/I;260PRINTUSING" #.###":Deviation(N,Sum,SumOfSquares)270NEXTI280END290!300EXTERNALSUBInit(N,Sum,SumOfSquares)310REM The parameters N, Sum, SumOfSquares are passed by reference. So, the arguments must be variables.320LETN=0330LETSum=0.0340LETSumOfSquares=0.0350ENDSUB360!370EXTERNALSUBAdd(N,Sum,SumOfSquares,P)380REM The parameters N, Sum, SumOfSquares are passed by reference. So, the arguments must be variables.390LETN=N+1400LETSum=Sum+P410LETSumOfSquares=SumOfSquares+P^2420ENDSUB430!440EXTERNALFUNCTIONDeviation(N,Sum,SumOfSquares)450LETDeviation=SQR(SumOfSquares/N-(Sum/N)^2)460ENDFUNCTION
N ITEM AVG STDDEV1 2 2.000 .0002 4 3.000 1.0003 4 3.333 .9434 4 3.500 .8665 5 3.800 .9806 5 4.000 1.0007 7 4.429 1.4008 9 5.000 2.000
Uses theMOD(array()) andSUM(array()) functions.
MAXITEMS=100FORi%=1TO8READnPRINT"Value = ";n", running SD = "FNrunningsd(n)NEXTENDDATA2,4,4,4,5,5,7,9DEFFNrunningsd(n)PRIVATElist(),i%DIMlist(MAXITEMS)i%+=1list(i%)=n=SQR(MOD(list())^2/i%-(SUM(list())/i%)^2)
Value = 2, running SD = 0Value = 4, running SD = 1Value = 4, running SD = 0.942809043Value = 4, running SD = 0.866025404Value = 5, running SD = 0.979795901Value = 5, running SD = 1Value = 7, running SD = 1.39970842Value = 9, running SD = 2
' FB 1.05.0 Win64FunctioncalcStandardDeviation(numberAsDouble)AsDoubleStatica()AsDoubleRedimPreservea(0ToUBound(a)+1)DimubAsUInteger=UBound(a)a(ub)=numberDimsumAsDouble=0.0ForiAsUInteger=0Toubsum+=a(i)NextDimmeanAsDouble=sum/(ub+1)DimdiffAsDoublesum=0.0ForiAsUInteger=0Toubdiff=a(i)-meansum+=diff*diffNextReturnSqr(sum/(ub+1))EndFunctionDima(0To7)AsDouble={2,4,4,4,5,5,7,9}ForiAsUInteger=0To7Print"Added";a(i);" SD now : ";calcStandardDeviation(a(i))NextPrintPrint"Press any key to quit"Sleep
Added 2 SD now : 0Added 4 SD now : 1Added 4 SD now : 0.9428090415820634Added 4 SD now : 0.8660254037844386Added 5 SD now : 0.9797958971132712Added 5 SD now : 1Added 7 SD now : 1.39970842444753Added 9 SD now : 2
100 PROGRAM "StDev.bas"110 LET N=8120 NUMERIC ARR(1 TO N)130 FOR I=1 TO N140 READ ARR(I)150 NEXT 160 DEF STDEV(N)170 LET S1,S2=0180 FOR I=1 TO N190 LET S1=S1+ARR(I)^2:LET S2=S2+ARR(I)200 NEXT 210 LET STDEV=SQR((N*S1-S2^2)/N^2)220 END DEF 230 FOR J=1 TO N240 PRINT J;"item =";ARR(J),"standard dev =";STDEV(J)250 NEXT 260 DATA 2,4,4,4,5,5,7,9
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 endfunction 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.000000New data 4 so S.D. now = 1.000000New data 4 so S.D. now = 0.942809New data 4 so S.D. now = 0.866025New data 5 so S.D. now = 0.979796New data 5 so S.D. now = 1.000000New data 7 so S.D. now = 1.399708New data 9 so S.D. now = 2.000000
;DefineourStandarddeviationfunctionDeclare.dStandard_deviation(x);MainprogramIfOpenConsole()Definei,xRestoreMyListFori=1To8Read.ixPrintN(StrD(Standard_deviation(x)))NextiPrint(#CRLF$+"Press ENTER to exit"):Input()EndIf;Calculationprocedure,withmemoryProcedure.dStandard_deviation(In)Staticin_summa,antalStaticin_kvadrater.qin_summa+inin_kvadrater+in*inantal+1ProcedureReturnPow((in_kvadrater/antal)-Pow(in_summa/antal,2),0.50)EndProcedure;datasectionDataSectionMyList:Data.i2,4,4,4,5,5,7,9EndDataSection
0.0000000000 1.0000000000 0.9428090416 0.8660254038 0.9797958971 1.0000000000 1.3997084244 2.0000000000
' Cumulative standard deviationTYPESampleNASINTEGERSumASSINGLESumOfSquaresASSINGLEENDTYPEDECLARESUBInit(SASSample)DECLARESUBAdd(SASSample,BYVALPASSINGLE)DECLAREFUNCTIONDeviation!(SASSample)DIMTest(1TO8)ASSINGLEDIMSASSampleFORI=1TO8READTest(I)NEXTIDATA2,4,4,4,5,5,7,9CALLInit(S)PRINT"N ITEM AVG STDDEV"FORI=LBOUND(Test)TOUBOUND(Test)CALLAdd(S,Test(I))PRINTUSING"# ";I;PRINTUSING" # ";Test(I);PRINTUSING"#.### ";S.Sum/I;PRINTUSING" #.###";Deviation!(S)NEXTIENDSUBAdd(SASSample,BYVALPASSINGLE)S.N=S.N+1S.Sum=S.Sum+PS.SumOfSquares=S.SumOfSquares+P^2ENDSUBFUNCTIONDeviation!(SASSample)Deviation!=SQR(S.SumOfSquares/S.N-(S.Sum/S.N)^2)ENDFUNCTIONSUBInit(SASSample)S.N=0S.Sum=0!S.SumOfSquares=0!ENDSUB
N ITEM AVG STDDEV1 2 2.000 0.0002 4 3.000 1.0003 4 3.333 0.9434 4 3.500 0.8665 5 3.800 0.9806 5 4.000 1.0007 7 4.429 1.4008 9 5.000 2.000
dim sdSave$(100) 'can call up to 100 versions 'holds (space-separated) number of data , sum of values and sum of squaressd$ = "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 num1 value in = 2 Stand Dev = 0.0000002 value in = 4 Stand Dev = 1.0000003 value in = 4 Stand Dev = 0.9428094 value in = 4 Stand Dev = 0.8660255 value in = 5 Stand Dev = 0.9797966 value in = 5 Stand Dev = 1.0000007 value in = 7 Stand Dev = 1.3997088 value in = 9 Stand Dev = 2.000000
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#include<stdio.h>#include<stdlib.h>#include<math.h>typedefenumAction{STDDEV,MEAN,VAR,COUNT}Action;typedefstructstat_obj_struct{doublesum,sum2;size_tnum;Actionaction;}sStatObject,*StatObject;StatObjectNewStatObject(Actionaction){StatObjectso;so=malloc(sizeof(sStatObject));so->sum=0.0;so->sum2=0.0;so->num=0;so->action=action;returnso;}#define FREE_STAT_OBJECT(so) \ free(so); so = NULLdoublestat_obj_value(StatObjectso,Actionaction){doublenum,mean,var,stddev;if(so->num==0.0)return0.0;num=so->num;if(action==COUNT)returnnum;mean=so->sum/num;if(action==MEAN)returnmean;var=so->sum2/num-mean*mean;if(action==VAR)returnvar;stddev=sqrt(var);if(action==STDDEV)returnstddev;return0;}doublestat_object_add(StatObjectso,doublev){so->num++;so->sum+=v;so->sum2+=v*v;returnstat_obj_value(so,so->action);}
doublev[]={2,4,4,4,5,5,7,9};intmain(){inti;StatObjectso=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);return0;}
usingSystem;usingSystem.Collections.Generic;usingSystem.Linq;namespacestandardDeviation{classProgram{staticvoidMain(string[]args){List<double>nums=newList<double>{2,4,4,4,5,5,7,9};for(inti=1;i<=nums.Count;i++)Console.WriteLine(sdev(nums.GetRange(0,i)));}staticdoublesdev(List<double>nums){List<double>store=newList<double>();foreach(doubleninnums)store.Add((n-nums.Average())*(n-nums.Average()));returnMath.Sqrt(store.Sum()/store.Count);}}}
010,9428090415820630,8660254037844390,97979589711327111,399708424447532
No attempt to handle different types -- standard deviation is intrinsically a real number.
#include<cassert>#include<cmath>#include<vector>#include<iostream>template<intN>structMomentsAccumulator_{std::vector<double>m_;MomentsAccumulator_():m_(N+1,0.0){}voidoperator()(doublev){doubleinc=1.0;for(auto&mi:m_){mi+=inc;inc*=v;}}};doubleStdev(conststd::vector<double>&moments){assert(moments.size()>2);assert(moments[0]>0.0);constdoublemean=moments[1]/moments[0];constdoublemeanSquare=moments[2]/moments[0];returnsqrt(meanSquare-mean*mean);}intmain(void){std::vector<int>data({2,4,4,4,5,5,7,9});MomentsAccumulator_<2>accum;for(autod:data){accum(d);std::cout<<"Running stdev: "<<Stdev(accum.m_)<<"\n";}}
(defnstateful-std-deviation[x](letfn[(std-dev[x](let[v(deref(find-var(symbol(str*ns*"/v"))))](swap!vconjx)(let[m(/(reduce+@v)(count@v))](Math/sqrt(/(reduce+(map#(*(-m%)(-m%))@v))(count@v))))))](when(nil?(resolve'v))(intern*ns*'v(atom[])))(std-devx)))
IDENTIFICATIONDIVISION.PROGRAM-ID.run-stddev.environmentdivision.input-outputsection.file-control. selectinput-fileassignto"input.txt" organizationislinesequential.data division.file section.fd input-file. 01inp-record. 03inp-fldpic 9(03).working-storagesection.01 fillerpic 9(01)value0. 88 no-more-inputvalue1.01 ws-tb-data. 03ws-tb-sizepic 9(03). 03ws-tb-table. 05ws-tb-fldpic s9(05)v9999comp-3occurs0to100timesdependingonws-tb-size.01 ws-stddevpic s9(05)v9999comp-3.PROCEDUREDIVISION. move0tows-tb-size openinputinput-file readinput-file atendsetno-more-inputtotrueend-readperformtestafteruntilno-more-inputadd1tows-tb-sizemoveinp-fldtows-tb-fld(ws-tb-size)call'stddev'usingbyreferencews-tb-dataws-stddevdisplay'inp='inp-fld' stddev='ws-stddevreadinput-fileatendsetno-more-inputtotrueend-readend-performcloseinput-file stoprun.end programrun-stddev.IDENTIFICATIONDIVISION.PROGRAM-ID.stddev.data division.working-storagesection.01 ws-tbxpic s9(03)comp.01 ws-tb-work. 03ws-sumpic s9(05)v9999comp-3value+0. 03ws-sumsqpic s9(05)v9999comp-3value+0. 03ws-avgpic s9(05)v9999comp-3value+0.linkagesection.01 ws-tb-data. 03ws-tb-sizepic 9(03). 03ws-tb-table. 05ws-tb-fldpic s9(05)v9999comp-3occurs0to100timesdependingonws-tb-size.01 ws-stddevpic s9(05)v9999comp-3.PROCEDUREDIVISIONusingws-tb-dataws-stddev. computews-sum=0performtestbeforevaryingws-tbxfrom1by+1untilws-tbx>ws-tb-sizecomputews-sum=ws-sum+ws-tb-fld(ws-tbx) end-perform computews-avgrounded=ws-sum/ws-tb-size computews-sumsq=0performtestbeforevaryingws-tbxfrom1by+1untilws-tbx>ws-tb-sizecomputews-sumsq=ws-sumsq+(ws-tb-fld(ws-tbx)-ws-avg)**2.0end-performcomputews-stddev=(ws-sumsq/ws-tb-size)**0.5goback.end programstddev.
sampleoutput:inp=002stddev=+00000.0000inp=004stddev=+00001.0000inp=004stddev=+00000.9427inp=004stddev=+00000.8660inp=005stddev=+00000.9797inp=005stddev=+00001.0000inp=007stddev=+00001.3996inp=009stddev=+00002.0000
Uses a class instance to maintain state.
classStandardDeviationconstructor:->@sum=0@sumOfSquares=0@values=0@deviation=0include:( n ) ->@values+=1@sum+=n@sumOfSquares+=n*nmean=@sum/@valuesmean*=mean@deviation=Math.sqrt@sumOfSquares/@values-meandev=newStandardDeviationvalues=[2,4,4,4,5,5,7,9]tmp=[]forvalueinvaluestmp.pushvaluedev.includevalueconsole.log""" Values:#{tmp} Standard deviation:#{dev.deviation} """
Values: 2Standard deviation: 0Values: 2,4Standard deviation: 1Values: 2,4,4Standard deviation: 0.9428090415820626Values: 2,4,4,4Standard deviation: 0.8660254037844386Values: 2,4,4,4,5Standard deviation: 0.9797958971132716Values: 2,4,4,4,5,5Standard deviation: 1Values: 2,4,4,4,5,5,7Standard deviation: 1.3997084244475297Values: 2,4,4,4,5,5,7,9Standard deviation: 2
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:
(defunrunning-stddev()(let((sum0)(sq0)(n0))(lambda(x)(incfsumx)(incfsq(*xx))(incfn)(/(sqrt(-(*nsq)(*sumsum)))n))))CL-USER>(loopwithf=(running-stddev)foriin'(24445579)do(formatt"~a ~a~%"i(funcallfi)))NIL20.041.040.9428090540.866025450.9797959351.071.399708592.0
In the REPL, one step at a time:
CL-USER>(setffn(running-stddev))#<InterpretedClosure(:INTERNALRUNNING-STDDEV)@#x21b9a492>CL-USER>(funcallfn2)0.0CL-USER>(funcallfn4)1.0CL-USER>(funcallfn4)0.94280905CL-USER>(funcallfn4)0.8660254CL-USER>(funcallfn5)0.97979593CL-USER>(funcallfn5)1.0CL-USER>(funcallfn7)1.3997085CL-USER>(funcallfn9)2.0
This example isincorrect. Please fix the code and remove this message.
|
BlackBox Component Builder
MODULEStandardDeviation;IMPORTStdLog,Args,Strings,Math;PROCEDUREMean(x:ARRAYOFREAL;n:INTEGER;OUTmean:REAL);VARi:INTEGER;total:REAL;BEGINtotal:=0.0;FORi:=0TOn-1DOtotal:=total+x[i]END;mean:=total/nENDMean;PROCEDURESDeviation(x:ARRAYOFREAL;n:INTEGER):REAL;VARi:INTEGER;mean,sum:REAL;BEGINMean(x,n,mean);sum:=0.0;FORi:=0TOn-1DOsum:=sum+((x[i]-mean)*(x[i]-mean));END;RETURNMath.Sqrt(sum/n);ENDSDeviation;PROCEDUREDo*;VARp:Args.Params;x:POINTERTOARRAYOFREAL;i,done:INTEGER;BEGINArgs.Get(p);IFp.argc>0THENNEW(x,p.argc);FORi:=0TOp.argc-1DOx[i]:=0.0END;FORi:=0TOp.argc-1DOStrings.StringToReal(p.args[i],x[i],done);StdLog.Int(i+1);StdLog.String(" :> ");StdLog.Real(SDeviation(x,i+1));StdLog.LnENDENDENDDo;ENDStandardDeviation.
Execute: ^Q StandardDeviation.Do 2 4 4 4 5 5 7 9 ~
1 :> 0.0 2 :> 1.0 3 :> 0.9428090415820634 4 :> 0.8660254037844386 5 :> 0.9797958971132712 6 :> 1.0 7 :> 1.39970842444753 8 :> 2.0
Use an object to keep state.
classStdDevAccumulatordefinitialize@n,@sum,@sum2=0,0.0,0.0enddef<<(num)@n+=1@sum+=num@sum2+=num**2Math.sqrt(@sum2*@n-@sum**2)/@n**2endendsd=StdDevAccumulator.newi=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.0adding 4: stddev of 2 samples is 1.0adding 4: stddev of 3 samples is 0.9428090415820634adding 4: stddev of 4 samples is 0.8660254037844386adding 5: stddev of 5 samples is 0.9797958971132712adding 5: stddev of 6 samples is 1.0adding 7: stddev of 7 samples is 1.3997084244475304adding 9: stddev of 8 samples is 2.0
defsdaccumn,sum,sum2=0,0.0,0.0->(num:Int32)don+=1sum+=numsum2+=num**2Math.sqrt((sum2*n-sum**2)/n**2)endendsd=sdaccum[2,4,4,4,5,5,7,9].each{|n|printsd.call(n),", "}
0.0, 1.0, 0.9428090415820634, 0.8660254037844386, 0.9797958971132712, 1.0, 1.3997084244475304, 2.0
importstd.stdio,std.math;structStdDev{realsum=0.0,sqSum=0.0;longnvalues;voidaddNumber(inrealinput)purenothrow{nvalues++;sum+=input;sqSum+=input^^2;}realgetStdDev()constpurenothrow{if(nvalues==0)return0.0;immutablerealmean=sum/nvalues;returnsqrt(sqSum/nvalues-mean^^2);}}voidmain(){StdDevstdev;foreach(el;[2.0,4,4,4,5,5,7,9]){stdev.addNumber(el);writefln("%e",stdev.getStdDev());}}
0.000000e+001.000000e+009.428090e-018.660254e-019.797959e-011.000000e+001.399708e+002.000000e+00
See:#Pascal
This entry is premised on the existence of a table, t, containing thenumbers (column x) in the required order, as specified by a columnof increasing values. The standard deviations are thencomputed on a rolling basis using DuckDB's support forincrementally growing "windows".
#CreateatablewiththedatatogetherwithasequencenumberCREATEORREPLACETABLEtASSELECTgenerate_subscripts(l,1)asid,unnest(l)::DOUBLEasx,from(values([2,4,4,4,5,5,7,9]))tbl(l);#ComputetherollingstandarddeviationswithouttheBesseladjustment:SELECT*,stddev_pop(x)OVER(ORDERBYid)as'cumulative stddev'FROMt;
┌───────┬────────┬────────────────────┐│ id │ x │ cumulative stddev ││ int64 │ double │ double │├───────┼────────┼────────────────────┤│ 1 │ 2.0 │ 0.0 ││ 2 │ 4.0 │ 1.0 ││ 3 │ 4.0 │ 0.9428090415820634 ││ 4 │ 4.0 │ 0.8660254037844386 ││ 5 │ 5.0 │ 0.9797958971132713 ││ 6 │ 5.0 │ 1.0 ││ 7 │ 7.0 │ 1.3997084244475304 ││ 8 │ 9.0 │ 2.0 │└───────┴────────┴────────────────────┘
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
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.01.00.94280904158206260.86602540378443860.97979589711327161.01.39970842444752972.0global sum sum2 n .func sd x . sum += x sum2 += x * x n += 1 return sqrt (sum2 / n - sum * sum / n / n).v[] = [ 2 4 4 4 5 5 7 9 ]for v in v[] print v & " " & sd v.
2 04 14 0.944 0.875 0.985 17 1.409 2
defmoduleStandard_deviationdodefadd_sample(pid,n),do:send(pid,{:add,n})defcreate,do:spawn_link(fn->loop([])end)defdestroy(pid),do:send(pid,:stop)defget(pid)dosend(pid,{:get,self()})receivedo{:get,value,_pid}->valueendenddeftaskdopid=create()forx<-[2,4,4,4,5,5,7,9],do:add_print(pid,x,add_sample(pid,x))destroy(pid)enddefpadd_print(pid,n,_add)doIO.puts"Standard deviation#{get(pid)} when adding#{n}"enddefploop(ns)doreceivedo{:add,n}->loop([n|ns]){:get,pid}->send(pid,{:get,loop_calculate(ns),self()})loop(ns):stop->:okendenddefploop_calculate(ns)doaverage=loop_calculate_average(ns):math.sqrt(loop_calculate_average(forx<-ns,do::math.pow(x-average,2)))enddefploop_calculate_average(ns),do:Enum.sum(ns)/length(ns)endStandard_deviation.task
Standard deviation 0.0 when adding 2Standard deviation 1.0 when adding 4Standard deviation 0.9428090415820634 when adding 4Standard deviation 0.8660254037844386 when adding 4Standard deviation 0.9797958971132712 when adding 5Standard deviation 1.0 when adding 5Standard deviation 1.3997084244475302 when adding 7Standard deviation 2.0 when adding 9
(defunrunning-std(items)(let((running-sum0)(running-len0)(running-squared-sum0)(result0))(dolist(itemitems)(setqrunning-sum(+running-sumitem))(setqrunning-len(1+running-len))(setqrunning-squared-sum(+running-squared-sum(*itemitem)))(setqresult(sqrt(-(/running-squared-sum(floatrunning-len))(/(*running-sumrunning-sum)(float(*running-lenrunning-len))))))(message"%f"result))result))(running-std'(24445579))
0.0000001.0000000.9428090.8660250.9797961.0000001.3997082.0000002.0
(let((x'(24445579)))(string-to-number(calc-eval"sqrt(vpvar($1))"nil(append'(vec)x))))
;; lexical-binding: t(require'generator)(iter-defunstd-dev-gen(lst)(let((sum0)(avg0)(tmp'())(std0))(dolist(ilst)(setqi(floati))(pushitmp)(setqsum(+sumi))(setqavg(/sum(lengthtmp)))(setqstd0)(dolist(jtmp)(setqstd(+std(expt(-javg)2))))(setqstd(/std(lengthtmp)))(setqstd(sqrtstd))(iter-yieldstd))))(let*((test-data'(24445579))(generator(std-dev-gentest-data)))(dolist(itest-data)(message"with %d: %f"i(iter-nextgenerator))))
-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}->Valueend.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->okend.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).
9> standard_deviation:task().Standard deviation 0.0 when adding 2Standard deviation 1.0 when adding 4Standard deviation 0.9428090415820634 when adding 4Standard deviation 0.8660254037844386 when adding 4Standard deviation 0.9797958971132712 when adding 5Standard deviation 1.0 when adding 5Standard deviation 1.3997084244475302 when adding 7Standard deviation 2.0 when adding 9
USING:accessorsiokernelmathmath.functionsmath.parsersequences;IN:standard-deviatorTUPLE:standard-deviatorsumsum^2n;:<standard-deviator>(--standard-deviator)0.0 0.0 0standard-deviatorboa;:current-std(standard-deviator--std)[[sum^2>>][n>>]bi/][[sum>>][n>>]bi/sq]bi-sqrt;:add-value(valuestandard-deviator--)[nip[1+]change-ndrop][[+]change-sumdrop][[[sq]dip+]change-sum^2drop]2tri;:main(--){2 4 4 4 5 5 7 9}<standard-deviator>[[add-value]curryeach]keepcurrent-stdnumber>stringprint;
01.01 C-- TEST SET01.10 S T(1)=2;S T(2)=4;S T(3)=4;S T(4)=401.20 S T(5)=5;S T(6)=5;S T(7)=7;S T(8)=901.30 D 2.101.35 T %6.4001.40 F I=1,8;S A=T(I);D 2.2;T "VAL",A;D 2.3;T " SD",A,!01.50 Q02.01 C-- RUNNING STDDEV02.02 C-- 2.1: INITIALIZE02.03 C-- 2.2: INSERT VALUE A02.04 C-- 2.3: A = CURRENT STDDEV02.10 S XN=0;S XS=0;S XQ=002.20 S XN=XN+1;S XS=XS+A;S XQ=XQ+A*A02.30 S A=FSQT(XQ/XN - (XS/XN)^2)
VAL= 2.00000 SD= 0.00000VAL= 4.00000 SD= 1.00000VAL= 4.00000 SD= 0.94281VAL= 4.00000 SD= 0.86603VAL= 5.00000 SD= 0.97980VAL= 5.00000 SD= 1.00000VAL= 7.00000 SD= 1.39971VAL= 9.00000 SD= 2.00000
:f+!( x addr -- )dupf@f+f!;:st-count( stats -- n )f@;:st-sum( stats -- sum )float+f@;:st-sumsq( stats -- sum*sum )2floats+f@;:st-mean( stats -- mean )dupst-sumst-countf/;:st-variance( stats -- var )dupst-sumsqdupst-meanfdupf*dupst-countf*f-st-countf/;:st-stddev( stats -- stddev )st-variancefsqrt;:st-add( fnum stats -- )dup1edupf+!float+fdupdupf+!float+fdupf*f+!std-stddev;
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 )2floats+f@;:st-variance( stats -- var )dupst-nvarst-countf/;:st-stddev( stats -- stddev )st-variancefsqrt;:st-add( x stats -- )dup1edupf+!\ update countfdupdupst-meanf-fswap( delta x )foverdupst-countf/( delta x delta/n )float+dupf+!\ update mean( delta x )dupf@f-f*float+f+!\ update nvarst-stddev;
Usage example:
createstats0ef,0ef,0ef,2estatsst-addf.\ 0.4estatsst-addf.\ 1.4estatsst-addf.\ 0.9428090415820634estatsst-addf.\ 0.8660254037844395estatsst-addf.\ 0.9797958971132715estatsst-addf.\ 1.7estatsst-addf.\ 1.399708424447539estatsst-addf.\ 2.
programstandard_deviationimplicit noneinteger(kind=4),parameter::dp=kind(0.0d0)real(kind=dp),dimension(:),allocatable::valsinteger(kind=4)::ireal(kind=dp),dimension(8)::sample_data=(/2,4,4,4,5,5,7,9/)doi=lbound(sample_data,1),ubound(sample_data,1)callsample_add(vals,sample_data(i))write(*,fmt='(''#'',I1,1X,''value = '',F3.1,1X,''stddev ='',1X,F10.8)')&i,sample_data(i),stddev(vals)end do if(allocated(vals))deallocate(vals)contains! Adds value :val: to array :population: dynamically resizing arraysubroutinesample_add(population,val)real(kind=dp),dimension(:),allocatable,intent(inout)::populationreal(kind=dp),intent(in)::valreal(kind=dp),dimension(:),allocatable::tmpinteger(kind=4)::nif(.not.allocated(population))then allocate(population(1))population(1)=valelsen=size(population)callmove_alloc(population,tmp)allocate(population(n+1))population(1:n)=tmppopulation(n+1)=valendif end subroutinesample_add! Calculates standard deviation for given set of valuesreal(kind=dp)functionstddev(vals)real(kind=dp),dimension(:),intent(in)::valsreal(kind=dp)::meaninteger(kind=4)::nn=size(vals)mean=sum(vals)/nstddev=sqrt(sum((vals-mean)**2)/n)end functionstddevend programstandard_deviation
#1 value = 2.0 stddev = 0.00000000#2 value = 4.0 stddev = 1.00000000#3 value = 4.0 stddev = 0.94280904#4 value = 4.0 stddev = 0.86602540#5 value = 5.0 stddev = 0.97979590#6 value = 5.0 stddev = 1.00000000#7 value = 7.0 stddev = 1.39970842#8 value = 9.0 stddev = 2.00000000
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 when specifying such values 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 being added in to V are never negative. This is demonstrated by comparing the results computed from StdDev(A), StdDev(A + 10), StdDev(A + 100), StdDev(A + 1000), etc.
Incidentally, Fortran implementations rarely enable re-entrancy 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.
REALFUNCTIONSTDDEV(X)!Standard deviation for successive values.REALX!The latest value.REALV!Scratchpad.INTEGERN!Ongoing: count of the values.REALEX,EX2!Ongoing: sum of X and X**2.SAVEN,EX,EX2!Retain values from one invocation to the next.DATAN,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.V=EX2/N-(EX/N)**2!The variance, but, it might come out negative!STDDEV=SIGN(SQRT(ABS(V)),V)!Protect the SQRT, but produce a negative result if so.END FUNCTIONSTDDEV!For the sequence of received X values.REALFUNCTIONSTDDEVP(X)!Standard deviation for successive values.REALX!The latest value.INTEGERN!Ongoing: count of the values.REALA,V!Ongoing: average, and sum of squared deviations.SAVEN,A,V!Retain values from one invocation to the next.DATAN,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 FUNCTIONSTDDEVP!For the sequence of received X values.REALFUNCTIONSTDDEVW(X)!Standard deviation for successive values.REALX!The latest value.REALV,D!Scratchpads.INTEGERN!Ongoing: count of the values.REALEX,EX2!Ongoing: sum of X and X**2.REALW!Ongoing: working mean.SAVEN,EX,EX2,W!Retain values from one invocation to the next.DATAN,EX,EX2/0,0.0,0.0/!Initial values.IF(N.LE.0)W=X!Take the first value as the working mean.N=N+1!Another value arrives.D=X-W!Its deviation from the working mean.EX=D+EX!Augment the total.EX2=D**2+EX2!Augment the sum of squares.V=EX2/N-(EX/N)**2!The variance, but, it might come out negative!STDDEVW=SIGN(SQRT(ABS(V)),V)!Protect the SQRT, but produce a negative result if so.END FUNCTIONSTDDEVW!For the sequence of received X values.REALFUNCTIONSTDDEVPW(X)!Standard deviation for successive values.REALX!The latest value.INTEGERN!Ongoing: count of the values.REALA,V!Ongoing: average, and sum of squared deviations.REALW!Ongoing: working mean.SAVEN,A,V,W!Retain values from one invocation to the next.DATAN,A,V/0,0.0,0.0/!Initial values.IF(N.LE.0)W=X!Oh for self-modifying code!N=N+1!Another value arrives.D=X-W!Its deviation from the working mean.V=(N-1)*(D-A)**2/N+V!First, as it requires the existing average.A=(D-A)/N+A!= [x + (n - 1).A)]/n: recover the total from the average.STDDEVPW=SQRT(V/N)!V can never be negative, even with limited precision.END FUNCTIONSTDDEVPW!For the sequence of received X values.PROGRAMTESTINTEGERI!A stepper.REALA(8)!The example data.DATAA/2.0,3*4.0,2*5.0,7.0,9.0/!Alas, another opportunity to use @ passed over.REALB!An offsetting base.WRITE(6,1)1FORMAT("Progressive calculation of the standard deviation."/1" I",7X,"A(I) EX EX2 Av V*N Ed Ed2 wAv V*N")B=1000000!Provoke truncation error.DOI=1,8!Step along the data series,WRITE(6,2)I,INT(A(I)+B),!No fractional part, so I don't want F11.0.1STDDEV(A(I)+B),STDDEVP(A(I)+B),!Showing progressive values.2STDDEVW(A(I)+B),STDDEVPW(A(I)+B)!These with a working mean.2FORMAT(I2,I11,1X,4F12.6)!Should do for the example.END DO!On to the next value.END
Output: the second pair of columns have the calculations done with a working mean and thus accumulate deviations from that.
Progressive calculation of the standard deviation.I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 2 0.000000 0.000000 0.000000 0.0000002 4 1.000000 1.000000 1.000000 1.0000003 4 0.942809 0.942809 0.942809 0.9428094 4 0.866025 0.866025 0.866025 0.8660255 5 0.979796 0.979796 0.979796 0.9797966 5 1.000000 1.000000 1.000000 1.0000007 7 1.399708 1.399708 1.399708 1.3997088 9 2.000000 2.000000 2.000000 2.000000
I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 12 0.000000 0.000000 0.000000 0.0000002 14 1.000000 1.000000 1.000000 1.0000003 14 0.942809 0.942809 0.942809 0.9428094 14 0.866025 0.866025 0.866025 0.8660255 15 0.979796 0.979796 0.979796 0.9797966 15 1.000000 1.000000 1.000000 1.0000007 17 1.399708 1.399708 1.399708 1.3997088 19 2.000000 2.000000 2.000000 2.000000
I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 102 0.000000 0.000000 0.000000 0.0000002 104 1.000000 1.000000 1.000000 1.0000003 104 0.942809 0.942809 0.942809 0.9428094 104 0.866025 0.866025 0.866025 0.8660255 105 0.979796 0.979796 0.979796 0.9797966 105 1.000000 0.999999 1.000000 1.0000007 107 1.399708 1.399708 1.399708 1.3997088 109 2.000000 1.999999 2.000000 2.000000
I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 1002 0.000000 0.000000 0.000000 0.0000002 1004 1.000000 1.000000 1.000000 1.0000003 1004 0.942809 0.942809 0.942809 0.9428094 1004 0.866025 0.866028 0.866025 0.8660255 1005 0.979796 0.979798 0.979796 0.9797966 1005 1.000000 1.000004 1.000000 1.0000007 1007 1.399708 1.399711 1.399708 1.3997088 1009 2.000000 1.999997 2.000000 2.000000
I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 10002 -2.000000 0.000000 0.000000 0.0000002 10004 -1.000000 1.000000 1.000000 1.0000003 10004 -0.666667 0.942809 0.942809 0.9428094 10004 1.936492 0.866072 0.866025 0.8660255 10005 2.181742 0.979829 0.979796 0.9797966 10005 2.309401 1.000060 1.000000 1.0000007 10007 1.801360 1.399745 1.399708 1.3997088 10009 2.645751 1.999987 2.000000 2.000000
I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 100002 19.493589 0.000000 0.000000 0.0000002 100004 7.416198 1.000000 1.000000 1.0000003 100004 -7.333333 0.942809 0.942809 0.9428094 100004 20.093531 0.865650 0.866025 0.8660255 100005 -1.280625 0.979531 0.979796 0.9797966 100005 -16.492422 1.000305 1.000000 1.0000007 100007 17.851427 1.399895 1.399708 1.3997088 100009 20.566963 1.999835 2.000000 2.000000
I A(I) EX EX2 Av V*N Ed Ed2 wAv V*N1 1000002 -80.024994 0.000000 0.000000 0.0000002 1000004 158.767120 1.000000 1.000000 1.0000003 1000004 -89.146576 0.942809 0.942809 0.9428094 1000004 90.795097 0.869074 0.866025 0.8660255 1000005 193.357590 0.981953 0.979796 0.9797966 1000005 238.361069 0.999691 1.000000 1.0000007 1000007 153.462296 1.399519 1.399708 1.3997088 1000009 151.284500 1.997653 2.000000 2.000000
Speaking loosely, to square a number of d digits accurately requires the ability to represent 2d digits accurately, with more digits needed if many such squares are to be added together accurately. In this example, 1000 when squared, is pushing at the limits of single precision. The average&variance method is resistant to this problem (and does not generate negative variances either!) because it manipulates differences from the running average, but it is still better to use a working mean.
In other words, a two-pass method will be more accurate (where the second pass calculates the variance by accumulating deviations from the actual average, itself calculated with a working mean) but at the cost of that second pass and the saving of all the values. Higher precision variables for the accumulations are the easiest way towards accurate results.
double local fn CalcSD( x as double ) static double n, sum, sum2 n++ sum += x sum2 += x * xend fn = sqr( sum2 / n - sum * sum / n / n )void local fn DoIt double testData(7) = {2,4,4,4,5,5,7,9} for int i = 0 to 7 double x = testData(i) double a = fn CalcSD( x ) printf @"value %.0f SD = %f", x, a nextend fnfn DoItHandleEventsvalue 2 SD = 0.000000value 4 SD = 1.000000value 4 SD = 0.942809value 4 SD = 0.866025value 5 SD = 0.979796value 5 SD = 1.000000value 7 SD = 1.399708value 9 SD = 2.000000
Algorithm to reduce rounding errors from WP article.
State maintained with a closure.
packagemainimport("fmt""math")funcnewRsdv()func(float64)float64{varn,a,qfloat64returnfunc(xfloat64)float64{n++a1:=a+(x-a)/nq,a=q+(x-a)*(x-a1),a1returnmath.Sqrt(q/n)}}funcmain(){r:=newRsdv()for_,x:=range[]float64{2,4,4,4,5,5,7,9}{fmt.Println(r(x))}}
010.94280904158206340.86602540378443860.979795897113271311.39970842444753042
Solution:
Listsamples=[]defstdDev={defsample->samples<<sampledefsum=samples.sum()defsumSq=samples.sum{it*it}defcount=samples.size()(sumSq/count - (sum/count)**2)**0.5}[2,4,4,4,5,5,7,9].each{println"${stdDev(it)}"}
010.94280904169991450.86602540378443860.979795897113271211.39970842434692622
We store the state in theST monad using anSTRef.
{-# LANGUAGE BangPatterns #-}importData.List(foldl')-- 'importData.STRefimportControl.Monad.STdataPairab=Pair!a!bsumLen::[Double]->PairDoubleDoublesumLen=fiof2.foldl'(\(Pairsl)x->Pair(s+x)(l+1))(Pair0.00)--'wherefiof2(Pairsl)=Pairs(fromIntegrall)divl::PairDoubleDouble->Doubledivl(Pair_0.0)=0.0divl(Pairsl)=s/lsd::[Double]->Doublesdxs=sqrt$foldl'(\ax->a+(x-m)^2)0xs/l--'wherep@(Pairsl)=sumLenxsm=divlpmkSD::STs(Double->STsDouble)mkSD=go<$>newSTRef[]wheregoaccx=modifySTRefacc(x:)>>(sd<$>readSTRefacc)main=mapM_print$runST$mkSD>>=forM[2.0,4.0,4.0,4.0,5.0,5.0,7.0,9.0]
Or, perhaps more simply, as a map-accumulation over an indexed list:
importData.List(mapAccumL)-------------- CUMULATIVE STANDARD DEVIATION -------------cumulativeStdDevns::[Float]->[Float]cumulativeStdDevns=snd.mapAccumLgo(0,0).zip[1.0..]wherego(s,q)(i,x)=let_s=s+x_q=q+(x^2)in((_s,_q),sqrt((_q/i)-((_s/i)^2)))--------------------------- TEST -------------------------main::IO()main=mapM_print$cumulativeStdDevns[2,4,4,4,5,5,7,9]
0.01.00.94280930.86602540.979795931.01.39970872.0
usingLambda;classMain{staticfunctionmain():Void{varnums=[2,4,4,4,5,5,7,9];for(iin1...nums.length+1)Sys.println(sdev(nums.slice(0,i)));}staticfunctionaverage<T:Float>(nums:Array<T>):Float{returnnums.fold(function(n,t)returnn+t,0)/nums.length;}staticfunctionsdev<T:Float>(nums:Array<T>):Float{varstore=[];varavg=average(nums);for(ninnums){store.push((n-avg)*(n-avg));}returnMath.sqrt(average(store));}}
010.9428090415820630.8660254037844390.97979589711327111.399708424447532
REAL :: n=8, set(n), sum=0, sum2=0set = (2,4,4,4,5,5,7,9)DO k = 1, n WRITE() 'Adding ' // set(k) // 'stdev = ' // stdev(set(k))ENDDOEND ! 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 = 0Adding 4 stdev = 1Adding 4 stdev = 0.9428090416Adding 4 stdev = 0.8660254038Adding 5 stdev = 0.9797958971Adding 5 stdev = 1Adding 7 stdev = 1.399708424Adding 9 stdev = 2
proceduremain()stddev()# reset state / emptyeverys:=stddev(![2,4,4,4,5,5,7,9])dowrite("stddev (so far) := ",s)endprocedurestddev(x)# running standard deviationstaticX,sumX,sum2Xif/xthen{# reset stateX:=[]sumX:=sum2X:=0.}else{# accumulateput(X,x)sumX+:=xsum2X+:=x^2returnsqrt((sum2X/*X)-(sumX/*X)^2)}end
stddev (so far) := 0.0stddev (so far) := 1.0stddev (so far) := 0.9428090415820626stddev (so far) := 0.8660254037844386stddev (so far) := 0.9797958971132716stddev (so far) := 1.0stddev (so far) := 1.39970842444753stddev (so far) := 2.0
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; seeMoving Average for an example. We will not reprise that here.
mean=:+/%#dev=:-meanstddevP=:[:%:@mean*:@devNB. A) 3 equivalent defs for stddevPstddevP=:[:mean&.:*:devNB. B) uses Under (&.:) to apply inverse of *: after meanstddevP=:%:@(mean@:*:-*:@mean)NB. C) sqrt of ((mean of squares) - (square of mean))stddevP\24445579010.9428090.8660250.97979611.399712
Alternatives:
Using verbose names for J primitives.
of=:@:sqrt=:%:sum=:+/squares=:*:data=:]mean=:sum%#stddevP=:sqrtofmeanofsquaresof(data-mean)stddevP\24445579010.9428090.8660250.97979611.399712
Or we could take a cue from theR implementation and reverse the Bessel correction to stddev:
require'stats'(%:@:(%~<:)@:#*stddev)\24445579010.9428090.8660250.97979611.399712
publicclassStdDev{intn=0;doublesum=0;doublesum2=0;publicdoublesd(doublex){n++;sum+=x;sum2+=x*x;returnMath.sqrt(sum2/n-sum*sum/n/n);}publicstaticvoidmain(String[]args){double[]testData={2,4,4,4,5,5,7,9};StdDevsd=newStdDev();for(doublex:testData){System.out.println(sd.sd(x));}}}
Uses a closure.
functionrunning_stddev(){varn=0;varsum=0.0;varsum_sq=0.0;returnfunction(num){n++;sum+=num;sum_sq+=num*num;returnMath.sqrt((sum_sq/n)-Math.pow(sum/n,2));}}varsd=running_stddev();varnums=[2,4,4,4,5,5,7,9];varstddev=[];for(variinnums)stddev.push(sd(nums[i]));// using WSHWScript.Echo(stddev.join(', ');
0, 1, 0.942809041582063, 0.866025403784439, 0.979795897113273, 1, 1.39970842444753, 2
Accumulating across a fold
(function(xs){returnxs.reduce(function(a,x,i){varn=i+1,sum_=a.sum+x,squaresSum_=a.squaresSum+(x*x);return{sum:sum_,squaresSum:squaresSum_,stages:a.stages.concat(Math.sqrt((squaresSum_/n)-Math.pow((sum_/n),2)))};},{sum:0,squaresSum:0,stages:[]}).stages})([2,4,4,4,5,5,7,9]);
[0,1,0.9428090415820626,0.8660254037844386,0.9797958971132716,1,1.3997084244475297,2]
As a map-accumulation:
(()=>{'use strict';// ---------- CUMULATIVE STANDARD DEVIATION ----------// cumulativeStdDevns :: [Float] -> [Float]constcumulativeStdDevns=ns=>{constgo=([s,q])=>([i,x])=>{const_s=s+x,_q=q+(x*x),j=1+i;return[[_s,_q],Math.sqrt((_q/j)-Math.pow(_s/j,2))];};returnmapAccumL(go)([0,0])(ns)[1];};// ---------------------- TEST -----------------------constmain=()=>showLog(cumulativeStdDevns([2,4,4,4,5,5,7,9]));// --------------------- GENERIC ---------------------// mapAccumL :: (acc -> x -> (acc, y)) -> acc -> [x] -> (acc, [y])constmapAccumL=f=>// A tuple of an accumulation and a list// obtained by a combined map and fold,// with accumulation from left to right.acc=>xs=>[...xs].reduce((a,x,i)=>{constpair=f(a[0])([i,x]);return[pair[0],a[1].concat(pair[1])];},[acc,[]]);// showLog :: a -> IO ()constshowLog=(...args)=>console.log(args.map(x=>JSON.stringify(x,null,2)).join(' -> '));// MAIN ---returnmain();})();
[ 0, 1, 0.9428090415820626, 0.8660254037844386, 0.9797958971132716, 1, 1.3997084244475297, 2]
We first define a filter, "simulate", that, if given a file ofobservations, will emit the standard deviations of the observationsseen 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)^2where 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:simulateExample 1
# observations.txt24445579
$jq-s-fDynamic_standard_deviation.jqobservations.txt010.94280904158206340.86602540378443860.97979589711327110.99999999999999991.39970842444753021.9999999999999998
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.4def simulate(stream): foreach stream as $observation (initial_state; update_state($observation); standard_deviation);
Example 2:
simulate( range(0;10) )
00.50.8164965809277261.1180339887498951.41421356237309511.70782512765993322.291287847477922.5819888974716112.8722813232690143
The following illustrates how jq can be used to process observations from an external (potentially unbounded) stream, one at a time. Here we usebash 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 PATHPROGRAM='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=nullwhileread-p"Next observation: "observationdoresult=$(echo"[$observation,$state ]"|jq-c"$PROGRAM")sed-n1p<<<"$result"state=$(sed-n2p<<<"$result")done
Example 3
$./standard_deviation_server.shNextobservation:100Nextobservation:205Nextobservation:08.16496580927726
Use a closure to create a running standard deviation function.
functionmakerunningstd(::Type{T}=Float64)whereT∑x=∑x²=zero(T)n=0functionrunningstd(x)∑x+=x∑x²+=x^2n+=1s=∑x²/n-(∑x/n)^2returnsendreturnrunningstdendtest=Float64[2,4,4,4,5,5,7,9]rstd=makerunningstd()println("Perform a running standard deviation of ",test)foriintestprintln(" - add$i → ",rstd(i))end
Perform a running standard deviation of [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0] - add 2.0 → 0.0 - add 4.0 → 1.0 - add 4.0 → 0.8888888888888875 - add 4.0 → 0.75 - add 5.0 → 0.9600000000000009 - add 5.0 → 1.0 - add 7.0 → 1.9591836734693864 - add 9.0 → 4.0
Using a class to keep the running sum, sum of squares and number of elements added so far:
// version 1.0.5-2classCumStdDev{privatevarn=0privatevarsum=0.0privatevarsum2=0.0funsd(x:Double):Double{n++sum+=xsum2+=x*xreturnMath.sqrt(sum2/n-sum*sum/n/n)}}funmain(args:Array<String>){valtestData=doubleArrayOf(2.0,4.0,4.0,4.0,5.0,5.0,7.0,9.0)valcsd=CumStdDev()for(dintestData)println("Add $d => ${csd.sd(d)}")}
Add 2.0 => 0.0Add 4.0 => 1.0Add 4.0 => 0.9428090415820626Add 4.0 => 0.8660254037844386Add 5.0 => 0.9797958971132708Add 5.0 => 1.0Add 7.0 => 1.399708424447531Add 9.0 => 2.0
// Stats computes a running mean and variance // See Knuth TAOCP vol 2, 3rd edition, page 232class Stats: M = 0.0 S = 0.0 n = 0 def incl(x): n += 1 if n == 1: M = x else: let mm = (x - M) M += mm / n S += mm * (x - M) def mean(): return M //def variance(): return (if n > 1.0: S / (n - 1.0) else: 0.0) // Bessel's correction def variance(): return (if n > 0.0: S / n else: 0.0) def stddev(): return sqrt(variance()) def count(): return n def test_stdv() -> float: let v = [2,4,4,4,5,5,7,9] let s = Stats {} for(v) x: s.incl(x+0.0) print concat_string(["Mean: ", string(s.mean()), ", Std.Deviation: ", string(s.stddev())], "")test_stdv()Mean: 5.0, Std.Deviation: 2.0
Uses a closure. Translation of JavaScript.
functionstdev()localsum,sumsq,k=0,0,0returnfunction(n)sum,sumsq,k=sum+n,sumsq+n^2,k+1returnmath.sqrt((sumsq/k)-(sum/k)^2)endendldev=stdev()fori,vinipairs{2,4,4,4,5,5,7,9}doprint(ldev(v))end
runningSTDDev[n_]:=(If[Not[ValueQ[$Data]],$Data={}];StandardDeviation[AppendTo[$Data,n]])
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 meanx2=cumsum(x.^2)./[1:n];% running squaresdev=sqrt(x2-m.*m)dev=0.000001.000000.942810.866030.979801.000001.399712.00000
Here is a vectorized one line solution as a function
functionstdDevEval(n)disp(sqrt(sum((n-sum(n)/length(n)).^2)/length(n)));end
StdDeviator={}StdDeviator.count=0StdDeviator.sum=0StdDeviator.sumOfSquares=0StdDeviator.add=function(x)self.count=self.count+1self.sum=self.sum+xself.sumOfSquares=self.sumOfSquares+x*xendfunctionStdDeviator.stddev=function()m=self.sum/self.countreturnsqrt(self.sumOfSquares/self.count-m*m)endfunctionsd=newStdDeviatorforxin[2,4,4,4,5,5,7,9]sd.addxendforprintsd.stddev
2
0П4П5П6С/ПП0ИП5+П5ИП0x^2ИП6+П6КИП4ИП6ИП4/ИП5ИП4/x^2-КвКорБП04
Instruction: В/О С/Пnumber С/Пnumber С/П ...
class StdDevdeclare ndeclare sumdeclare sum2def StdDev()n = 0sum = 0sum2 = 0enddef sd(x)this.n += 1this.sum += xthis.sum2 += x*xreturn sqrt(sum2/n - sum*sum/n/n)endendtestData = {2,4,4,4,5,5,7,9}sd = new(StdDev)for x in testDataprintln sd.sd(x)end0.01.00.94280904158206340.86602540378443860.97979589711327121.01.39970842444753042.0
importmath,strutilsvarsdSum,sdSum2,sdN=0.0procsd(x:float):float=sdN+=1sdSum+=xsdSum2+=x*xsqrt(sdSum2/sdN-sdSum*sdSum/(sdN*sdN))forvaluein[float2,4,4,4,5,5,7,9]:echovalue," ",formatFloat(sd(value),precision=-1)
2 04 14 0.9428094 0.8660255 0.9797965 17 1.399719 2
importmath,strutilstypeSDAccum=objectsdN,sdSum,sdSum2:floatvaraccum:SDAccumprocadd(accum:varSDAccum;value:float):float=# Add a value to the accumulator. Return the standard deviation.accum.sdN+=1accum.sdSum+=valueaccum.sdSum2+=value*valueresult=sqrt(accum.sdSum2/accum.sdN-accum.sdSum*accum.sdSum/(accum.sdN*accum.sdN))forvaluein[float2,4,4,4,5,5,7,9]:echovalue," ",formatFloat(accum.add(value),precision=-1)
Same output.
importmath,strutilsfuncaccumBuilder():auto=varsdSum,sdSum2,sdN=0.0result=func(value:float):float=sdN+=1sdSum+=valuesdSum2+=value*valueresult=sqrt(sdSum2/sdN-sdSum*sdSum/(sdN*sdN))letstd=accumBuilder()forvaluein[float2,4,4,4,5,5,7,9]:echovalue," ",formatFloat(std(value),precision=-1)
Same output.
MODULECumulativeStandardDeviation;IMPORTMath,Out;TYPECSD*=POINTERTOCSDStruct;CSDStruct=RECORDsum,sum2,sd:REAL;n:INTEGEREND;PROCEDUREsample*(sd:CSD;x:REAL):REAL;VARr,sum,sum2:REAL;BEGINsum:=sd.sum+x;sum2:=sd.sum2+(x*x);;INC(sd.n);r:=FLT(sd.n);IFsd.n=0THENsd.sd:=0.0ELSEsd.sd:=Math.sqrt(sum2/r-sum*sum/r/r)END;sd.sum:=sum;sd.sum2:=sum2RETURNsd.sdENDsample;PROCEDURENewCSD*():CSD;VARsd:CSD;BEGINNEW(sd);sd.n:=0;sd.sum:=0.0;sd.sum2:=0.0;sd.sd:=0.0RETURNsdENDNewCSD;PROCEDUREshowSampleEffect(sd:CSD;x:REAL);VARdev:REAL;BEGINdev:=sample(sd,x);Out.Int(sd.n,6);Out.String(": ");Out.Real(x,10);Out.String(" -> ");Out.Real(dev,10);Out.LnENDshowSampleEffect;PROCEDUREtest;VARsd:CSD;BEGINsd:=NewCSD();showSampleEffect(sd,2.0);showSampleEffect(sd,4.0);showSampleEffect(sd,4.0);showSampleEffect(sd,4.0);showSampleEffect(sd,5.0);showSampleEffect(sd,5.0);showSampleEffect(sd,7.0);showSampleEffect(sd,9.0)ENDtest;BEGINtest;ENDCumulativeStandardDeviation.
Showing the output fromoberonc, the real number format may be different with other compilers.
1: 2.000000 -> 0.000000 2: 4.000000 -> 1.000000 3: 4.000000 -> 0.942809 4: 4.000000 -> 0.866025 5: 5.000000 -> 0.979796 6: 5.000000 -> 1.000000 7: 7.000000 -> 1.399708 8: 9.000000 -> 2.000000
Oberon-2 has type bound procedures, allowing us to writesd.sample( x ) instead ofsample( sd, x ).
MODULECumulativeStandardDeviation;IMPORTMath,Out;TYPECSD*=POINTERTOCSDStruct;CSDStruct*=RECORDsum,sum2,sd:REAL;n:INTEGER;END;PROCEDURE(sd:CSD)sample*(x:REAL):REAL;VARr,sum,sum2:REAL;BEGINsum:=sd.sum+x;sum2:=sd.sum2+(x*x);;INC(sd.n);r:=sd.n;IFsd.n=0THENsd.sd:=0.0ELSEsd.sd:=Math.sqrt(sum2/r-sum*sum/r/r)END;sd.sum:=sum;sd.sum2:=sum2;RETURNsd.sdENDsample;PROCEDURENewCSD*():CSD;VARsd:CSD;BEGINNEW(sd);sd.n:=0;sd.sum:=0.0;sd.sum2:=0.0;sd.sd:=0.0;RETURNsdENDNewCSD;PROCEDURE(sd:CSD)showSampleEffect(x:REAL);VARdev:REAL;BEGINdev:=sd.sample(x);Out.Int(sd.n,6);Out.String(": ");Out.Real(x,15);Out.String(" -> ");Out.Real(dev,15);Out.LnENDshowSampleEffect;PROCEDUREtest;VARsd:CSD;BEGINsd:=NewCSD();sd.showSampleEffect(2.0);sd.showSampleEffect(4.0);sd.showSampleEffect(4.0);sd.showSampleEffect(4.0);sd.showSampleEffect(5.0);sd.showSampleEffect(5.0);sd.showSampleEffect(7.0);sd.showSampleEffect(9.0)ENDtest;BEGINtest;ENDCumulativeStandardDeviation.
Similar to theOberon-07 sample, except for differences in real number formatting.
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(); } }}#import <Foundation/Foundation.h>@interfaceSDAccum :NSObject{doublesum,sum2;unsignedintnum;}-(double)value:(double)v;-(unsignedint)count;-(double)mean;-(double)variance;-(double)stddev;@end@implementationSDAccum-(double)value:(double)v{sum+=v;sum2+=v*v;num++;return[selfstddev];}-(unsignedint)count{returnnum;}-(double)mean{return(num>0)?sum/(double)num:0.0;}-(double)variance{doublem=[selfmean];return(num>0)?(sum2/(double)num-m*m):0.0;}-(double)stddev{returnsqrt([selfvariance]);}@endintmain(){@autoreleasepool{doublev[]={2,4,4,4,5,5,7,9};SDAccum*sdacc=[[SDAccumalloc]init];for(inti=0;i<sizeof(v)/sizeof(*v);i++)printf("adding %f\tstddev = %f\n",v[i],[sdaccvalue:v[i]]);}return0;}
#import <Foundation/Foundation.h>typedefdouble(^Func)(double);// a block that takes a double and returns a doubleFuncsdCreator(){__blockintn=0;__blockdoublesum=0;__blockdoublesum2=0;return^(doublex){sum+=x;sum2+=x*x;n++;returnsqrt(sum2/n-sum*sum/n/n);};}intmain(){@autoreleasepool{doublev[]={2,4,4,4,5,5,7,9};Funcsdacc=sdCreator();for(inti=0;i<sizeof(v)/sizeof(*v);i++)printf("adding %f\tstddev = %f\n",v[i],sdacc(v[i]));}return0;}
letsqrx=x*.xletstddevl=letn,sx,sx2=List.fold_left(fun(n,sx,sx2)x->succn,sx+.x,sx2+.sqrx)(0,0.,0.)linsqrt((sx2-.sqrsx/.floatn)/.floatn)let_=letl=[2.;4.;4.;4.;5.;5.;7.;9.]inPrintf.printf"List: ";List.iter(Printf.printf"%g ")l;Printf.printf"\nStandard deviation: %g\n"(stddevl)
List: 2 4 4 4 5 5 7 9Standard deviation: 2
Oforth does not have global variables that can be used to create statefull functions.
Here, we create a channel to hold current list of numbers. Constraint is that this channel can't hold mutable objects. On the other hand, stddev function is thread safe and can be called by tasks running in parallel.
Channel new [ ] over send drop const: StdValues: stddev(x)| l | StdValues receive x + dup ->l StdValues send drop #qs l map sum l size asFloat / l avg sq - sqrt ;
>[ 2, 4, 4, 4, 5, 5, 7, 9 ] apply(#[ stddev println ])010.9428090415820630.8660254037844390.97979589711327211.399708424447532ok>
sdacc=.SDAccum~newx=.array~of(2,4,4,4,5,5,7,9)sd=0doi=1tox~sizesd=sdacc~value(x[i])Say'#'i'value ='x[i]'stdev ='sdend::classSDAccum::methodsumattribute::methodsum2attribute::methodcountattribute::methodinitself~sum=0.0self~sum2=0.0self~count=0::methodvalueexposesumsum2countparseargxsum=sum+xsum2=sum2+x*xcount=count+1returnself~stddev::methodmeanexposesumcountreturnsum/count::methodvarianceexposesum2countm=self~meanreturnsum2/count-m*m::methodstddevreturnself~sqrt(self~variance)::methodsqrtargnifn=0thenreturn0ans=n/2prev=ndountilprev=ansprev=ansans=(prev+(n/prev))/2endreturnans
#1 value = 2 stdev = 0#2 value = 4 stdev = 1#3 value = 4 stdev = 0.94280905#4 value = 4 stdev = 0.866025405#5 value = 5 stdev = 0.979795895#6 value = 5 stdev = 1#7 value = 7 stdev = 1.39970844#8 value = 9 stdev = 2
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])programstddev;usesmath;constn=8;vararr:array[1..n]ofreal=(2,4,4,4,5,5,7,9);functionstddev(n:integer):real;vari:integer;s1,s2,variance,x:real;beginfori:=1tondobeginx:=arr[i];s1:=s1+power(x,2);s2:=s2+xend;variance:=((n*s1)-(power(s2,2)))/(power(n,2));stddev:=sqrt(variance)end;vari:integer;beginfori:=1tondobeginwriteln(i,' item=',arr[i]:2:0,' stddev=',stddev(i):18:15)endend.
1 item= 2 stddev= 0.0000000000000002 item= 4 stddev= 1.0000000000000003 item= 4 stddev= 0.9428090415820644 item= 4 stddev= 0.8660254037844395 item= 5 stddev= 0.9797958971132716 item= 5 stddev= 1.0000000000000007 item= 7 stddev= 1.3997084244475308 item= 9 stddev= 2.000000000000000
programprj_CalcStdDerv;{$APPTYPE CONSOLE}usesMath;varSeries:ArrayofExtended;UserString:String;functionAppendAndCalc(NewVal:Extended):Extended;beginsetlength(Series,high(Series)+2);Series[high(Series)]:=NewVal;result:=PopnStdDev(Series);end;constdata:array[0..7]ofExtended=(2,4,4,4,5,5,7,9);varrr:Extended;beginsetlength(Series,0);forrrindatadobeginwriteln(rr,' -> ',AppendAndCalc(rr));end;Readln;end.
2.0000000000000000E+0000 -> 0.0000000000000000E+0000 4.0000000000000000E+0000 -> 1.0000000000000000E+0000 4.0000000000000000E+0000 -> 9.4280904158206337E-0001 4.0000000000000000E+0000 -> 8.6602540378443865E-0001 5.0000000000000000E+0000 -> 9.7979589711327124E-0001 5.0000000000000000E+0000 -> 1.0000000000000000E+0000 7.0000000000000000E+0000 -> 1.3997084244475303E+0000 9.0000000000000000E+0000 -> 2.0000000000000000E+0000
##functionsdcreator():real->real;beginvarSum:=0.0;varSum2:=0.0;varN:=0;result:=function(x:real):real->beginN+=1;Sum+=x;Sum2+=x*x;result:=sqrt(Sum2/N-Sum*Sum/(N*N));end;end;varsd:=sdcreator();foreachvarvaluein|2,4,4,4,5,5,7,9|doprintln(value,sd(value));
2 04 14 0.9428090415820644 0.8660254037844395 0.9797958971132725 17 1.399708424447539 2
{packageSDAccum;subnew{my$class=shift;my$self={};$self->{sum}=0.0;$self->{sum2}=0.0;$self->{num}=0;bless$self,$class;return$self;}subcount{my$self=shift;return$self->{num};}submean{my$self=shift;return($self->{num}>0)?$self->{sum}/$self->{num}:0.0;}subvariance{my$self=shift;my$m=$self->mean;return($self->{num}>0)?$self->{sum2}/$self->{num}-$m*$m:0.0;}substddev{my$self=shift;returnsqrt($self->variance);}subvalue{my$self=shift;my$v=shift;$self->{sum}+=$v;$self->{sum2}+=$v*$v;$self->{num}++;return$self->stddev;}}
my$sdacc=SDAccum->new;my$sd;foreachmy$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;substddev{my$x=shift;$num++;returnsqrt(($sum2+=$x**2)/$num-(($sum+=$x)/$num)**2);}}printstddev($_),"\n"forqw(2 4 4 4 5 5 7 9);
010.9428090415820630.8660254037844390.97979589711327211.399708424447532
one-liner:
perl-MMath::StdDev-e'$d=new Math::StdDev;foreach my $v ( 2,4,4,4,5,5,7,9 ) {$d->Update($v); print $d->variance(),"\n"}'small script:
useMath::StdDev;$d=newMath::StdDev;foreachmy$v(2,4,4,4,5,5,7,9){$d->Update($v);print$d->variance(),"\n"}
010.9428090415820630.8660254037844390.97979589711327111.399708424447532
demo\rosetta\Standard_deviation.exw contains a copy of this code and a version that could be the basis for a library version that can handle multiple active data sets concurrently.
withjavascript_semanticsatomsdn=0,sdsum=0,sdsumsq=0proceduresdadd(atomn)sdn+=1sdsum+=nsdsumsq+=n*nendprocedurefunctionsdavg()returnsdsum/sdnendfunctionfunctionsddev()returnsqrt(sdsumsq/sdn-power(sdsum/sdn,2))endfunction--test code:constanttestset={2,4,4,4,5,5,7,9}integertifori=1tolength(testset)doti=testset[i]sdadd(ti)printf(1,"N=%d Item=%d Avg=%5.3f StdDev=%5.3f\n",{i,ti,sdavg(),sddev()})endfor
N=1 Item=2 Avg=2.000 StdDev=0.000N=2 Item=4 Avg=3.000 StdDev=1.000N=3 Item=4 Avg=3.333 StdDev=0.943N=4 Item=4 Avg=3.500 StdDev=0.866N=5 Item=5 Avg=3.800 StdDev=0.980N=6 Item=5 Avg=4.000 StdDev=1.000N=7 Item=7 Avg=4.429 StdDev=1.400N=8 Item=9 Avg=5.000 StdDev=2.000
As of 1.0.6 there is a builtin std_dev(), the following produces the same output as above, albeit w/o any running totals:
withjavascript_semanticsconstanttestset={2,4,4,4,5,5,7,9}fori,tintestsetdosequencetest=testset[1..i]printf(1,"N=%d Item=%d Avg=%5.3f StdDev=%5.3f\n",{i,t,average(test),std_dev(test)})endfor
This is just straight PHP class usage, respecting the specifications "stateful" and "one at a time":
<?phpclasssdcalc{private$cnt,$sumup,$square;function__construct(){$this->reset();}# callable on an instancefunctionreset(){$this->cnt=0;$this->sumup=0;$this->square=0;}functionadd($f){$this->cnt++;$this->sumup+=$f;$this->square+=pow($f,2);return$this->calc();}functioncalc(){if($this->cnt==0||$this->sumup==0){return0;}else{returnsqrt($this->square/$this->cnt-pow(($this->sumup/$this->cnt),2));}}}# start test, adding test data one by one$c=newsdcalc();foreach([2,4,4,4,5,5,7,9]as$v){printf('Adding %g: result %g%s',$v,$c->add($v),PHP_EOL);}
This will produce the output:
Adding 2: result 0Adding 4: result 1Adding 4: result 0.942809Adding 4: result 0.866025Adding 5: result 0.979796Adding 5: result 1Adding 7: result 1.39971Adding 9: result 2
(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)) ) )
2.00 -> 0.004.00 -> 1.004.00 -> 0.944.00 -> 0.875.00 -> 0.985.00 -> 1.007.00 -> 1.409.00 -> 2.00
*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;AVERAGE= 5.50000E+0000;Standard deviation 2.87228E+0000
require"table2"localfmt=require"fmt"localfunctioncum_std_dev(a)fori=1,#adolocalb=a:slice(1,i)fmt.print("Values : %,s",b)coroutine.yield(b:stddev())endendlocala={2,4,4,4,5,5,7,9}localco=coroutine.create(cum_std_dev)whiletruedolocal_,sd=coroutine.resume(co,a)if!sdthenbreakendfmt.print("Std Dev : %10.8f\n",sd)end
Values : {2}Std Dev : 0.00000000Values : {2, 4}Std Dev : 1.00000000Values : {2, 4, 4}Std Dev : 0.94280904Values : {2, 4, 4, 4}Std Dev : 0.86602540Values : {2, 4, 4, 4, 5}Std Dev : 0.97979590Values : {2, 4, 4, 4, 5, 5}Std Dev : 1.00000000Values : {2, 4, 4, 4, 5, 5, 7}Std Dev : 1.39970842Values : {2, 4, 4, 4, 5, 5, 7, 9}Std Dev : 2.00000000This implementation takes the form of an advanced function which can act like a cmdlet and receive input from the pipeline.
functionGet-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-StandardDeviation010.9428090415820630.8660254037844390.97979589711327111.399708424447532
The program should work with Python 2.x and 3.x, although the output would not be a tuple in 3.x
>>>frommathimportsqrt>>>defsd(x):sd.sum+=xsd.sum2+=x*xsd.n+=1.0sum,sum2,n=sd.sum,sd.sum2,sd.nreturnsqrt(sum2/n-sum*sum/n/n)>>>sd.sum=sd.sum2=sd.n=0>>>forvaluein(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)>>>
>>>classSD(object):# Plain () for python 3.xdef__init__(self):self.sum,self.sum2,self.n=(0,0,0)defsd(self,x):self.sum+=xself.sum2+=x*xself.n+=1.0sum,sum2,n=self.sum,self.sum2,self.nreturnsqrt(sum2/n-sum*sum/n/n)>>>sd_inst=SD()>>>forvaluein(2,4,4,4,5,5,7,9):print(value,sd_inst.sd(value))
You could rename the methodsd to__call__ this would make the class instance callable like a function so instead of usingsd_inst.sd(value) it would change tosd_inst(value) for the same results.
>>>frommathimportsqrt>>>defsdcreator():sum=sum2=n=0defsd(x):nonlocalsum,sum2,nsum+=xsum2+=x*xn+=1.0returnsqrt(sum2/n-sum*sum/n/n)returnsd>>>sd=sdcreator()>>>forvaluein(2,4,4,4,5,5,7,9):print(value,sd(value))20.041.040.94280904158240.86602540378450.97979589711351.071.3997084244592.0
>>>frommathimportsqrt>>>defsdcreator():sum=sum2=n=0whileTrue:x=yieldsqrt(sum2/n-sum*sum/n/n)ifnelseNonesum+=xsum2+=x*xn+=1.0>>>sd=sdcreator()>>>sd.send(None)>>>forvaluein(2,4,4,4,5,5,7,9):print(value,sd.send(value))20.041.040.94280904158240.86602540378450.97979589711351.071.3997084244592.0
>>>myMean=lambdaMyList:reduce(lambdax,y:x+y,MyList)/float(len(MyList))>>>myStd=lambdaMyList:(reduce(lambdax,y:x+y,map(lambdax:(x-myMean(MyList))**2,MyList))/float(len(MyList)))**.5>>>printmyStd([2,4,4,4,5,5,7,9])2.0
To compute the running sum, one must keep track of the number of items, the sum of values, and the sum of squares.
If the goal is to get a vector of running standard deviations, the simplest is to do it with cumsum:
cumsd<-function(x){n<-seq_along(x)sqrt(cumsum(x^2)/n-(cumsum(x)/n)^2)}set.seed(12345L)x<-rnorm(10)cumsd(x)# [1] 0.0000000 0.3380816 0.8752973 1.1783628 1.2345538 1.3757142 1.2867220 1.2229056 1.1665168 1.1096814# Compare to the naive implementation, i.e. compute sd on each sublist:Vectorize(function(k)sd(x[1:k])*sqrt((k-1)/k))(seq_along(x))# [1] NA 0.3380816 0.8752973 1.1783628 1.2345538 1.3757142 1.2867220 1.2229056 1.1665168 1.1096814# Note that the first is NA because sd is unbiased formula, hence there is a division by n-1, which is 0 for n=1.
The task requires an accumulator solution:
accumsd<-function(){n<-0m<-0s<-0function(x){n<<-n+1m<<-m+xs<<-s+x*xsqrt(s/n-(m/n)^2)}}f<-accumsd()sapply(x,f)# [1] 0.0000000 0.3380816 0.8752973 1.1783628 1.2345538 1.3757142 1.2867220 1.2229056 1.1665168 1.1096814
#langracket(requiremath)(definerunning-stddev(let([ns'()])(λ(n)(set!ns(consnns))(stddevns))));; run it on each number, return the last result(last(maprunning-stddev'(24445579)))
(formerly Perl 6)
Using a closure:
subsd (@a) {my$mean =@aR/ [+]@a;sqrt@aR/ [+]map (* -$mean)²,@a;}subsdaccum {my@a;return {push@a,$^x;sd@a; };}my&f =sdaccum;sayf$_for2,4,4,4,5,5,7,9;
Using a state variable (remember that<(x-<x>)²> = <x²> - <x>²):
substddev($x) {sqrt ( .[2] +=$x²) / ++.[0] - ((.[1] +=$x ) / .[0])²givenstate @;}say .&stddevfor<2 4 4 4 5 5 7 9>;
010.9428090415820630.8660254037844390.97979589711327111.399708424447532
These REXX versions use running sums.
/*REXX program calculates and displays the standard deviation of a given set of numbers.*/parsearg#/*obtain optional arguments from the CL*/if#=''then#=24445579/*None specified? Then use the default*/n=words(#);$=0;$$=0;L=length(n)/*N: # items; $,$$: sums to be zeroed*//* [↓] process each number in the list*/doj=1forn_=word(#,j);$=$+_$$=$$+_**2say' item'right(j,L)":"right(_,4)' average='left($/j,12),' standard deviation='sqrt($$/j-($/j)**2)end/*j*//* [↑] prettify output with whitespace*/say'standard deviation: 'sqrt($$/n-($/n)**2)/*calculate & display the std deviation*/exit0/*stick a fork in it, we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/sqrt:procedure;parseargx;ifx=0thenreturn0;d=digits();h=d+6;m.=9;numericformnumericdigits;parsevalueformat(x,2,1,,0)'E0'withg'E'_.;g=g*.5'e'_%2doj=0whileh>9;m.j=h;h=h%2+1;end/*j*/dok=j+5to0by-1;numericdigitsm.k;g=(g+x/g)*.5;end/*k*/numericdigitsd;returng/1
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= 2standard deviation: 2
/*REXX program calculates and displays the standard deviation of a given set of numbers.*/parsearg#/*obtain optional arguments from the CL*/if#=''then#=24445579/*None specified? Then use the default*/n=words(#);$=0;$$=0/*N: # items; $,$$: sums to be zeroed*//* [↓] process each number in the list*/doj=1forn/*perform summation on two sets of #'s.*/_=word(#,j);$=$+_$$=$$+_**2end/*j*/say'standard deviation: 'sqrt($$/n-($/n)**2)/*calculate&display the std, deviation.*/exit0/*stick a fork in it, we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/sqrt:procedure;parseargx;ifx=0thenreturn0;d=digits();h=d+6;m.=9;numericformnumericdigits;parsevalueformat(x,2,1,,0)'E0'withg'E'_.;g=g*.5'e'_%2doj=0whileh>9;m.j=h;h=h%2+1;end/*j*/dok=j+5to0by-1;numericdigitsm.k;g=(g+x/g)*.5;end/*k*/numericdigitsd;returng/1
standard deviation: 2
# Project : Cumulative standard deviationdecimals(6)sdsave = list(100) sd = "2,4,4,4,5,5,7,9"sumval = 0sumsqs = 0for num = 1 to 8 sd = substr(sd, ",", "") stddata = number(sd[num]) sumval = sumval + stddata sumsqs = sumsqs + pow(stddata,2) standdev = pow(((sumsqs / num) - pow((sumval /num),2)),0.5) sdsave[num] = string(num) + " " + string(sumval) +" " + string(sumsqs) see "" + num + " value in = " + stddata + " Stand Dev = " + standdev + nlnext
Output:
1 value in = 2 Stand Dev = 02 value in = 4 Stand Dev = 13 value in = 4 Stand Dev = 0.9428094 value in = 4 Stand Dev = 0.8660255 value in = 5 Stand Dev = 0.9797966 value in = 5 Stand Dev = 17 value in = 7 Stand Dev = 1.3997088 value in = 9 Stand Dev = 2
≪ CL∑ { } SWAP 1 OVER SIZEFOR j DUP j GET ∑+IF j 1 >THEN SDEV ∑DAT SIZE 1 GET DUP 1 - SWAP / √ * ROT SWAP + SWAPENDNEXT DROP CL∑ ≫ 'CSDEV' STO≪ CL∑ 1 ≪ ∑+ PSDEV ≫ DOSUBS CL∑≫ 'CSDEV' STO1: { 0 1 0.942809041582 0.866025403784 0.979795897113 1 1.39970842445 2 }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.
classStdDevAccumulatordefinitialize@n,@sum,@sumofsquares=0,0.0,0.0enddef<<(num)# return self to make this possible: sd << 1 << 2 << 3 # => 0.816496580927726@n+=1@sum+=num@sumofsquares+=num**2selfenddefstddevMath.sqrt((@sumofsquares/@n)-(@sum/@n)**2)enddefto_sstddev.to_sendendsd=StdDevAccumulator.newi=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.0adding 4: stddev of 2 samples is 1.0adding 4: stddev of 3 samples is 0.942809041582063adding 4: stddev of 4 samples is 0.866025403784439adding 5: stddev of 5 samples is 0.979795897113272adding 5: stddev of 6 samples is 1.0adding 7: stddev of 7 samples is 1.39970842444753adding 9: stddev of 8 samples is 2.0
defsdaccumn,sum,sum2=0,0.0,0.0lambdado|num|n+=1sum+=numsum2+=num**2Math.sqrt((sum2/n)-(sum/n)**2)endendsd=sdaccum[2,4,4,4,5,5,7,9].each{|n|printsd.call(n),", "}
0.0, 1.0, 0.942809041582063, 0.866025403784439, 0.979795897113272, 1.0, 1.39970842444753, 2.0,
Using a struct:
pub struct CumulativeStandardDeviation { n: f64, sum: f64, sum_sq: f64}impl CumulativeStandardDeviation { pub fn new() -> Self { CumulativeStandardDeviation { n: 0., sum: 0., sum_sq: 0. } } fn push(&mut self, x: f64) -> f64 { self.n += 1.; self.sum += x; self.sum_sq += x * x; (self.sum_sq / self.n - self.sum * self.sum / self.n / self.n).sqrt() }}fn main() { let nums = [2, 4, 4, 4, 5, 5, 7, 9]; let mut cum_stdev = CumulativeStandardDeviation::new(); for num in nums.iter() { println!("{}", cum_stdev.push(*num as f64)); }}010.94280904158206260.86602540378443860.979795897113270811.3997084244475312
Using a closure:
fn sd_creator() -> impl FnMut(f64) -> f64 { let mut n = 0.0; let mut sum = 0.0; let mut sum_sq = 0.0; move |x| { sum += x; sum_sq += x*x; n += 1.0; (sum_sq / n - sum * sum / n / n).sqrt() }}fn main() { let nums = [2, 4, 4, 4, 5, 5, 7, 9]; let mut sd_acc = sd_creator(); for num in nums.iter() { println!("{}", sd_acc(*num as f64)); }}010.94280904158206260.86602540378443860.979795897113270811.3997084244475312
*--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;
N SD1 0.000002 1.000003 0.942814 0.866035 0.979806 1.000007 1.399718 2.00000
import scala.math.sqrtobject StddevCalc extends App { def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = { def avg(ts: Iterable[T])(implicit num: Fractional[T]): T = num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T val mean: T = avg(ts) // Leave val type of T // Root of mean diffs val stdDev = sqrt(ts.map { x => val diff = num.toDouble(num.minus(x, mean)) diff * diff }.sum / ts.size) (mean, stdDev) } 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)))) println(s"Successfully completed without errors. [total ${scala.compat.Platform.currentTime - executionStart}ms]")}(define (standart-deviation-generator) (let ((nums '())) (lambda (x) (set! nums (cons x nums)) (let* ((mean (/ (apply + nums) (length nums))) (mean-sqr (lambda (y) (expt (- y mean) 2))) (variance (/ (apply + (map mean-sqr nums)) (length nums)))) (sqrt variance)))))(let loop ((f (standart-deviation-generator)) (input '(2 4 4 4 5 5 7 9))) (unless (null? input) (display (f (car input))) (newline) (loop f (cdr input))))
Scilab has the built-in functionstdev 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.
T=[2,4,4,4,5,5,7,9];stdev(T)*sqrt((length(T)-1)/length(T))
-->T=[2,4,4,4,5,5,7,9];-->stdev(T)*sqrt((length(T)-1)/length(T)) ans = 2.
Using an object to keep state:
class StdDevAccumulator(n=0, sum=0, sumofsquares=0) { method <<(num) { n += 1 sum += num sumofsquares += num**2 self } method stddev { sqrt(sumofsquares/n - pow(sum/n, 2)) } method to_s { self.stddev.to_s }} var i = 0var sd = StdDevAccumulator()[2,4,4,4,5,5,7,9].each {|n| say "adding #{n}: stddev of #{i+=1} samples is #{sd << n}"}adding 2: stddev of 1 samples is 0adding 4: stddev of 2 samples is 1adding 4: stddev of 3 samples is 0.942809041582063365867792482806465385713114583585adding 4: stddev of 4 samples is 0.866025403784438646763723170752936183471402626905adding 5: stddev of 5 samples is 0.979795897113271239278913629882356556786378992263adding 5: stddev of 6 samples is 1adding 7: stddev of 7 samples is 1.39970842444753034182701947126050936683768427466adding 9: stddev of 8 samples is 2
Usingstatic variables:
func stddev(x) { static(num=0, sum=0, sum2=0) num++ sqrt( (sum2 += x**2) / num - (((sum += x) / num)**2) )} %n(2 4 4 4 5 5 7 9).each { say stddev(_) }010.9428090415820633658677924828064653857131145835850.8660254037844386467637231707529361834714026269050.97979589711327123927891362988235655678637899226311.399708424447530341827019471260509366837684274662
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.-- the minimal tablecreate table if not exists teststd (n double precision not null);-- code modularity with view, we could have used a common table expression insteadcreate view vteststd as select count(n) as cnt, sum(n) as tsum, sum(power(n,2)) as tsqrfrom teststd;-- you can of course put this code into every querycreate 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,9insert 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 datadelete from teststd;
With a command likepsql <rosetta-std-dev.sql you will get an output like this: (duplicate lines generously deleted, locale is DE)
CREATE TABLEFEHLER: Relation »vteststd« existiert bereitsCREATE FUNCTIONINSERT 0 1 std_deviation --------------- 0(1 Zeile)INSERT 0 1 std_deviation --------------- 1 0.942809041582063 0.866025403784439 0.979795897113272 1 1.39970842444753 2DELETE 8
import Darwinclass 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()value 2 SD = 0.0value 4 SD = 1.0value 4 SD = 0.942809041582063value 4 SD = 0.866025403784439value 5 SD = 0.979795897113271value 5 SD = 1.0value 7 SD = 1.39970842444753value 9 SD = 2.0
Functional:
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.8742514835862standardDeviation([2,4,4,4,5,5,7,9]) // 2.0or
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])} }}# Demonstrationset sdacc [SDAccum new]foreach val {2 4 4 4 5 5 7 9} { set sd [$sdacc value $val]}puts "the standard deviation is: $sd"the standard deviation is: 2.0
# Make a coroutine out of a lambda applicationcoroutine 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)}] }}}# Demonstrationforeach val {2 4 4 4 5 5 7 9} { set sd [sd $val]}sd stopputs "the standard deviation is: $sd"NOTE. The helper functionavg is not namedaverage to avoid a name conflict withWorksheetFunction.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 FunctionFunction 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 FunctionmakeArr: If 9 = Err.Number Then ReDim list(0) Else 'something's wrong Err.Raise Err.Number End If Resume NextEnd FunctionSub 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)) NextEnd Sub
0 1 0.942809041582063 0.866025403784439 0.979795897113271 1 1.39970842444753 2
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.WriteLineNextFunction sd(arr,n)mean = 0variance = 0For j = 0 To nmean = mean + arr(j)Nextmean = mean/(n+1)For k = 0 To nvariance = variance + ((arr(k)-mean)^2)Nextvariance = variance/(n+1)sd = FormatNumber(Sqr(variance),6)End Function
value = 2 running sd = 0.000000value = 4 running sd = 1.000000value = 4 running sd = 0.942809value = 4 running sd = 0.866025value = 5 running sd = 0.979796value = 5 running sd = 1.000000value = 7 running sd = 1.399708value = 9 running sd = 2.000000
import "./fmt" for Fmtimport "./math" for Numsvar cumStdDev = Fiber.new { |a| for (i in 0...a.count) { var b = a[0..i] System.print("Values : %(b)") Fiber.yield(Nums.popStdDev(b)) }}var a = [2, 4, 4, 4, 5, 5, 7, 9]while (true) { var sd = cumStdDev.call(a) if (cumStdDev.isDone) return Fmt.print("Std Dev : $10.8f\n", sd)}Values : [2]Std Dev : 0.00000000Values : [2, 4]Std Dev : 1.00000000Values : [2, 4, 4]Std Dev : 0.94280904Values : [2, 4, 4, 4]Std Dev : 0.86602540Values : [2, 4, 4, 4, 5]Std Dev : 0.97979590Values : [2, 4, 4, 4, 5, 5]Std Dev : 1.00000000Values : [2, 4, 4, 4, 5, 5, 7]Std Dev : 1.39970842Values : [2, 4, 4, 4, 5, 5, 7, 9]Std Dev : 2.00000000
include c:\cxpl\codes; \intrinsic 'code' declarationsint 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);]
0.00000 1.00000 0.94281 0.86603 0.97980 1.00000 1.39971 2.00000
Works with: 0.14.x, 0.15.0-dev.849+36499c251
const std = @import("std");pub fn main() !void { const stdout = std.io.getStdOut().writer(); const arr = [_]f32{ 2, 4, 4, 4, 5, 5, 7, 9 }; var sum: f32 = 0; var sum2: f32 = 0; var count: f32 = 0; for (arr) |n| { sum += n; sum2 += std.math.pow(f32, n, 2); count += 1; const cumulative: f32 = count * sum2 - std.math.pow(f32, sum, 2); const running_sd: f32 = std.math.sqrt(cumulative) / count; try stdout.print("{d} {d}\n", .{n, running_sd}); }}2 04 14 0.942809054 0.86602545 0.979795935 17 1.39970859 2
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())2zkl: sd:=sdf(); sd(2);sd(4);sd(4);sd(4);sd(5);sd(5);sd(7);sd(9)2