It's likely my old stuff was in APL.
But for now, I've worked up the following code which should work for
several
groups of data, stored in three columns containing numbers, means and
unbiased
standard deviations.
(Later thought while drafting - I think it should work for biased estimates
too, if we replace n - 1, below, by n - UNBIASED where global UNBIASED is
1 or 0 respectively. Recode accordingly if you wish!)
NB. retrieve sum of sqs of x given n, m, s for one set
sx2 =: 3 : 0
'n m s' =. y
(n - 1 0) (+/ . *) *: s, m
)
NB. sd given n, m, sx2
sd =: 3 : 0
'n m s2' =. y NB. number, mean, sum of squared values
%: (n - 1) %~ s2 - n (* *:) m
)
NB. Given n1, m1, s1, n2, m2, s2, etc, in 3 columns,
NB. what is overall n,m,s for two sets combined?
nmscomb =: 3 : 0
'n m s' =. |:y
sumxx =. +/sx2"1 y NB. overall sum of squares of values
nmall =. ({:, %/) n +/ . * m,. 1 NB. overall number & mean
(}:, sd) nmall, sumxx NB. append overall standard
deviation
:
nmscomb x,:y NB. for two sets in lh & rh args
)
5 51.4 23.6072 nmscomb 12 70.4167 26.3972 NB. Devon's "data" as 2 sets
17 64.8236 26.4226
2 3 $ 5 51.4 23.6072 12 70.4167 26.3972 NB. as table of with rows
5 51.4 23.6072
12 70.4167 26.3972
nmscomb 2 3 $ 5 51.4 23.6072 12 70.4167 26.3972 NB. process table
17 64.8236 26.4226
Any use?
Mike
On 20/09/2018 18:01, 'Mike Day' via Programming wrote:
I was going to have a look back at my old workings from when I was
actually
involved professionally in thinking about meta-stats, but you've
evidently
found something useful while we were having supper!
Let me know if you think it might still be useful, always assuming I
can find
anything relevant!
Cheers,
Mike
[NB no snipping of back-thread!]
On 20/09/2018 17:20, Devon McCormick wrote:
OK - it turns out that there's a very simple solution, courtesy of this
site:
https://www.mathworks.com/matlabcentral/fileexchange/35605-overall-mean-standard-deviation-of-groups-of-observations?w.mathworks.com
.
combineSDsScalars=: 3 : 0
'grpn grpmean grpstd'=. y
alln=. +/grpn NB. Total # observations
allmean=. grpmean +/ . * grpn%alln NB. Weighted mean of all
ESS=. (*:grpstd) +/ . * <:grpn NB. Total error sum of squares
GSS=. grpn +/ . * *:grpmean-allmean NB. Total group sum of squares
%:(ESS+GSS)%<:alln NB. Total standard deviation
)
nn=. 23 81 52 67 34; 97 72 53 22 85 39 97 94 82 35 79 90
<"1|:(#,mean,stddev)&>nn
+----+------------+---------------+
|5 12|51.4 70.4167|23.6072 26.3972|
+----+------------+---------------+
combineSDsScalars <"1|:(#,mean,stddev)&>nn
26.4226
stddev ;nn NB. Correct answer
26.4226
nn2=. 59 52 90 62 0; 64 29 14 30 92 98 82 12 54 30; 84 55 35 7 8
23 66
60 99 62 93 91 35 41 27
stddev ;nn2 NB. Correct answer
30.2466
combineSDsScalars <"1|:(#,mean,stddev)&>nn2
30.2466
On Thu, Sep 20, 2018 at 11:40 AM Raul Miller <[email protected]>
wrote:
I think the approach should extend to multiple series if you build a
combine routine which works in a reduce context.
In other words, you need to build all the summary statistics when
combining two series. Or, something like this (warning: untested
code):
combineSum=: +&{.
combineMean=: +/&(*/)&}: % combineSum
combine=: combineSum, combineMean, combineSDsS@,@,.
I *think* combineSDsS would be like combineSDsScalars, but without the
m argument (which doesn't seem to be relevant). If that's right, you'd
use it like this:
combine&;/summary1;summary2;summary3
I hope this helps,
--
Raul
On Thu, Sep 20, 2018 at 11:18 AM Devon McCormick <[email protected]>
wrote:
I work with a system that handles a lot of complicated problems very
simply
but I've recently run into a difficult problem, the solution of which
would
be how do we re-construct the overall standard deviation of a
number of
series, to which we would rather not refer to in detail, when we have
only
the sample size, the mean, and the standard deviation for each group.
How do we combine these?
There's a lot of stuff on the web but most of what I've found is
wrong -
unless I'm wrong and the numerous different answers are all right -
but
for
this reference - almost:
http://atozmath.com/CONM/Ch2_CombinedSD.aspx?q1=1.1%605%2c12%6051.4%2c70.4167%6023.6072%2c26.3972%60SD2#PrevPart
.
I say "almost" because, on the website, the formula is correct but the
arithmetic is wrong!
The site shows an example worked out in painstaking detail - given the
summary stats for this example - but the final answer is shown as
26.1952
when the correct answer is 26.4226.
However, if I properly reproduce the formula and do my own math, I get
the
correct answer, to whit - here's the formula, as given, with the math,
and
following is my own work in J, reproducing the formulas but getting
the
math right.
Combined Mean = 64.8236
Combined Standard deviation :
σ12=√(N1-1)⋅σ21+(N2-1)⋅σ22+N1⋅N2N1+N2⋅(ˉx21+ˉx22-2ˉx1ˉx2)N1+N2-1
]nn=. (5?@$100);12?@$100
+--------------+-----------------------------------+
|23 81 52 67 34|97 72 53 22 85 39 97 94 82 35 79 90|
+--------------+-----------------------------------+
stddev ;nn NB. Correct answer
26.4226
(*:@:stddev)&>nn
557.3 696.811
var&>nn
557.3 696.811
(#,mean,stddev)&>nn
5 51.4 23.6072
12 70.4167 26.3972
NB. Following based on formula here:
http://atozmath.com/CONM/Ch2_CombinedSD.aspx?q1=1.1%605%2c12%6051.4%2c70.4167%6023.6072%2c26.3972%60SD2#PrevPart
combineSDs=: [: %: ((([: <: #&>) +/ .* var&>) + ([: (*/ % +/)
#&>) *
([:
+/ [: *: mean&>) - [: +: [: */ mean&>) % [: <: [: +/ #&>
combineSDs nn NB. This works based on the mean,
variance, and
counts of nn.
26.4226
NB. This version makes it clear we can do this using only the
summary
statistics:
'n1 n2'=. #&>nn [ 's21 s22'=. var&>nn [ 'm1 m2'=. mean&>nn [
m=. mean
;nn
%:(((n1-1)*s21)+((n2-1)*s22)+((n1*n2)%n1+n2)*(*:m1)+(*:m2)-+:m1*m2)%<:n1+n2
26.4226
combineSDsScalars=: 3 : 0
'n1 n2 m1 m2 m s21 s22'=. y
%:(((n1-1)*s21)+((n2-1)*s22)+((n1*n2)%n1+n2)*(*:m1)+(*:m2)-+:m1*m2)%<:n1+n2
)
combineSDsScalars n1,n2,m1,m2,m,s21,s22
26.4226
But how do I extend this for more than two series?
--
Devon McCormick, CFA
Quantitative Consultant
----------------------------------------------------------------------
For information about J forums see
http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
