Statistics/Basic
Statistics is all about large groups of numbers. When talking about a set of sampled data, most frequently used is their mean value and standard deviation (stddev). If you have set of data where , the mean is , while the stddev is .
When examining a large quantity of data, one often uses a histogram, which shows the counts of data samples falling into a prechosen set of intervals (or bins). When plotted, often as bar graphs, it visually indicates how often each data value occurs.
Task Using your language's random number routine, generate real numbers in the range of [0, 1]. It doesn't matter if you chose to use open or closed range. Create 100 of such numbers (i.e. sample size 100) and calculate their mean and stddev. Do so for sample size of 1,000 and 10,000, maybe even higher if you feel like. Show a histogram of any of these sets. Do you notice some patterns about the standard deviation?
Extra Sometimes so much data need to be processed that it's impossible to keep all of them at once. Can you calculate the mean, stddev and histogram of a trillion numbers? (You don't really need to do a trillion numbers, just show how it can be done.)
- Hint
For a finite population with equal probabilities at all points, one can derive:
C
Sample code. <lang C>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
- include <stdint.h>
- define n_bins 10
double rand01() { return rand() / (RAND_MAX + 1.0); }
double avg(int count, double *stddev, int *hist) { double x[count]; double m = 0, s = 0;
for (int i = 0; i < n_bins; i++) hist[i] = 0; for (int i = 0; i < count; i++) { m += (x[i] = rand01()); hist[(int)(x[i] * n_bins)] ++; }
m /= count; for (int i = 0; i < count; i++) s += x[i] * x[i]; *stddev = sqrt(s / count - m * m); return m; }
void hist_plot(int *hist) { int max = 0, step = 1; double inc = 1.0 / n_bins;
for (int i = 0; i < n_bins; i++) if (hist[i] > max) max = hist[i];
/* scale if numbers are too big */ if (max >= 60) step = (max + 59) / 60;
for (int i = 0; i < n_bins; i++) { printf("[%5.2g,%5.2g]%5d ", i * inc, (i + 1) * inc, hist[i]); for (int j = 0; j < hist[i]; j += step) printf("#"); printf("\n"); } }
/* record for moving average and stddev. Values kept are sums and sum data^2
* to avoid excessive precision loss due to divisions, but some loss is inevitable */
typedef struct { uint64_t size; double sum, x2; uint64_t hist[n_bins]; } moving_rec;
void moving_avg(moving_rec *rec, double *data, int count) { double sum = 0, x2 = 0; /* not adding data directly to the sum in case both recorded sum and * count of this batch are large; slightly less likely to lose precision*/ for (int i = 0; i < count; i++) { sum += data[i]; x2 += data[i] * data[i]; rec->hist[(int)(data[i] * n_bins)]++; }
rec->sum += sum; rec->x2 += x2; rec->size += count; }
int main() { double m, stddev; int hist[n_bins], samples = 10;
while (samples <= 10000) { m = avg(samples, &stddev, hist); printf("size %5d: %g %g\n", samples, m, stddev); samples *= 10; }
printf("\nHistograph:\n"); hist_plot(hist);
printf("\nMoving average:\n N Mean Sigma\n"); moving_rec rec = { 0, 0, 0, {0} }; double data[100]; for (int i = 0; i < 10000; i++) { for (int j = 0; j < 100; j++) data[j] = rand01();
moving_avg(&rec, data, 100);
if ((i % 1000) == 999) { printf("%4lluk %f %f\n", rec.size/1000, rec.sum / rec.size, sqrt(rec.x2 * rec.size - rec.sum * rec.sum)/rec.size ); } } }</lang>
Go
<lang go>package main
import (
"fmt" "math" "rand" "strings"
)
func main() {
sample(100) sample(1000) sample(10000)
}
func sample(n int) {
// generate data
d := make([]float64, n)
for i := range d {
d[i] = rand.Float64()
}
// show mean, standard deviation
var sum, ssq float64
for _, s := range d {
sum += s
ssq += s * s
}
fmt.Println(n, "numbers")
m := sum / float64(n)
fmt.Println("Mean: ", m)
fmt.Println("Stddev:", math.Sqrt(ssq/float64(n)-m*m))
// show histogram
h := make([]int, 10)
for _, s := range d {
h[int(s*10)]++
}
for _, c := range h {
fmt.Println(strings.Repeat("*", c*205/int(n)))
}
fmt.Println()
}</lang> Output:
100 numbers Mean: 0.5231064889267764 Stddev: 0.292668237816841 **************** **************** ************************ ********************** ****************** ****************** **************** ************************** ************************ ******************** 1000 numbers Mean: 0.496026080160094 Stddev: 0.2880988956436907 ********************* ******************** ***************** *********************** ****************** ********************** ******************** ********************* ****************** ******************* 10000 numbers Mean: 0.5009091903581223 Stddev: 0.289269693719711 ******************* ******************** ******************** ******************** ********************* ******************** ******************* ******************* ******************** *********************
For the extra, here is an outline of a map reduce strategy. The main task indicated that numbers should be generated before doing any computations on them. Consistent with that, The function getSegment returns data based on a starting and ending index, as if it were accessing some large data store. <lang go>package main
import (
"fmt" "math" "rand" "strings"
)
func main() {
bigSample(1e7)
}
func bigSample(n int64) {
sum, ssq, h := reduce(0, n)
// compute final statistics and output as above
fmt.Println(n, "numbers")
m := sum / float64(n)
fmt.Println("Mean: ", m)
fmt.Println("Stddev:", math.Sqrt(ssq/float64(n)-m*m))
for _, c := range h {
fmt.Println(strings.Repeat("*", c*205/int(n)))
}
fmt.Println()
}
const threshold = 1e6
func reduce(start, end int64) (sum, ssq float64, h []int) {
n := end - start
if n < threshold {
d := getSegment(start, end)
return computeSegment(d)
}
// map to two sub problems
half := (start + end) / 2
sum1, ssq1, h1 := reduce(start, half)
sum2, ssq2, h2 := reduce(half, end)
// combine results
for i, c := range h2 {
h1[i] += c
}
return sum1 + sum2, ssq1 + ssq2, h1
}
func getSegment(start, end int64) []float64 {
d := make([]float64, end-start)
for i := range d {
d[i] = rand.Float64()
}
return d
}
func computeSegment(d []float64) (sum, ssq float64, h []int) {
for _, s := range d {
sum += s
ssq += s * s
}
h = make([]int, 10)
for _, s := range d {
h[int(s*10)]++
}
return
}</lang> Output:
10000000 numbers Mean: 0.4999673191148989 Stddev: 0.2886663876567514 ******************** ******************** ******************** ******************** ******************** ******************** ******************** ******************** ******************** ********************
Icon and Unicon
The following uses the stddev procedure from the Standard_deviation task. In this example,
<lang Icon>procedure main(A)
W := 50 # avg width for histogram bar B := 10 # histogram bins if *A = 0 then put(A,100) # 100 if none specified
while N := get(A) do { # once per argument
write("\nN=",N)
N := 0 < integer(N) | next # skip if invalid
stddev() # reset
m := 0.
H := list(B,0) # Histogram of
every i := 1 to N do { # calc running ...
s := stddev(r := ?0) # ... std dev
m +:= r/N # ... mean
H[integer(*H*r)+1] +:= 1 # ... histogram
}
write("mean=",m)
write("stddev=",s)
every i := 1 to *H do # show histogram
write(right(real(i)/*H,5)," : ",repl("*",integer(*H*50./N*H[i])))
}
end</lang>
Output:
N=100 mean=0.4941076275054806 stddev=0.2812938788216594 0.1 : **************************************** 0.2 : ******************************************************* 0.3 : ******************************************************* 0.4 : ********************************************************************** 0.5 : **************************************** 0.6 : ********************************************* 0.7 : **************************************** 0.8 : ***************************************************************** 0.9 : **************************************** 1.0 : ************************************************** N=10000 mean=0.4935428224375008 stddev=0.2884171825227816 0.1 : *************************************************** 0.2 : *************************************************** 0.3 : *************************************************** 0.4 : ************************************************** 0.5 : **************************************************** 0.6 : ************************************************* 0.7 : *********************************************** 0.8 : ************************************************ 0.9 : ************************************************** 1.0 : *********************************************** N=1000000 mean=0.4997503773607869 stddev=0.2886322440610256 0.1 : ************************************************* 0.2 : ************************************************** 0.3 : ************************************************** 0.4 : ************************************************** 0.5 : ************************************************* 0.6 : ************************************************** 0.7 : ************************************************* 0.8 : ************************************************* 0.9 : ************************************************** 1.0 : *************************************************
J
J has library routines to compute mean and standard deviation: <lang j> require'statfns'
(mean,stddev) ?1000#0
0.484669 0.287482
(mean,stddev) ?10000#0
0.503642 0.290777
(mean,stddev) ?100000#0
0.499677 0.288726</lang>
And, for a histogram:
<lang j>histogram=: <: @ (#/.~) @ (i.@#@[ , I.) require'plot' plot ((%*1+i.)100) ([;histogram) ?10000#0</lang>
but these are not quite what is being asked for here.
Instead:
<lang j>histogram=: <: @ (#/.~) @ (i.@#@[ , I.)
meanstddevP=:3 :0
NB. compute mean and std dev of y random numbers NB. picked from even distribution between 0 and 1 NB. and display a normalized ascii histogram for this sample NB. note: should use population mean, not sample mean, for stddev NB. given the equation specified for this task. h=.s=.t=. 0 buckets=. (%~1+i.)10 for_n.i.<.y%1e6 do. data=. ?1e6#0 h=.h+ buckets histogram data s=.s+ +/ data t=.t+ +/(data-0.5)^2 end. data=. ?(1e6|y)#0 h=.h+ buckets histogram data s=.s+ +/ data t=.t++/(data-0.5)^2 smoutput (<.300*h%y)#"0'#' (s%y),%:t%y
)</lang>
Example use:
<lang j> meanstddevP 1000
0.488441 0.289744
meanstddevP 10000
0.49697 0.289433
meanstddevP 100000
0.500872 0.288241</lang>
(That said, note that these numbers are random, so reported standard deviation will vary with the random sample being tested.)
This could handle a trillion random numbers on a bog-standard computer, but I am not inclined to wait that long.
Liberty BASIC
Be aware that the PRNG in LB has a SLIGHT bias. <lang lb> call sample 100 call sample 1000 call sample 10000
end
sub sample n
dim dat( n)
for i =1 to n
dat( i) =rnd( 1)
next i
'// show mean, standard deviation
sum =0
sSq =0
for i =1 to n
sum =sum +dat( i)
sSq =sSq +dat( i)^2
next i
print n; " data terms used."
mean =sum / n print "Mean ="; mean
print "Stddev ="; ( sSq /n -mean^2)^0.5
'// show histogram
nBins =10
dim bins( nBins)
for i =1 to n
z =int( nBins *dat( i))
bins( z) =bins( z) +1
next i
for b =0 to nBins -1
for j =1 to int( nBins *bins( b)) /n *70)
print "#";
next j
print
next b
print
end sub </lang>
100000 data terms used. Mean =0.49870232 Stddev =0.28926563 ###################################################################### ###################################################################### ###################################################################### ###################################################################### ##################################################################### ##################################################################### ##################################################################### ##################################################################### ###################################################################### #####################################################################
PARI/GP
<lang parigp>mean(v)={
sum(i=1,#v,v[i])/#v
}; stdev(v,mu="")={
if(mu=="",mu=mean(v)); sqrt(sum(i=1,#v,(v[i]-mu)^2))/#v
}; histogram(v,bins=16,low=0,high=1)={
my(u=vector(bins),width=(high-low)/bins); for(i=1,#v,u[(v[i]-low)\width+1]++); u
}; show(n)={
my(v=vector(n,i,random(1.)),mu=mean(v),s=stdev(v,mu),h=histogram(v),sz=ceil(n/50/16));
for(i=1,16,for(j=1,h[i]\sz,print1("#"));print());
print("Mean: "mu);
print("Stdev: "s);
}; show(100); show(1000); show(10000);</lang>
Perl
<lang perl>my @histogram = (); my $hist_rows = 10; my $sum = 0; my $sum_squares = 0; my $n = $ARGV[ 0 ];
for( my $counter = 0; $counter < $n; $counter ++ ){
my $current = rand( 1 ); $sum += $current; $sum_squares += $current ** 2; $histogram[ $current * $hist_rows ] += 1;
}
my $mean = $sum / $n;
print "$n numbers\n" .
"Mean: $mean\n" .
"Stddev: " . sqrt(($sum_squares / $n) - ($mean ** 2)) . "\n";
for( my $row_counter = 0; $row_counter < $hist_rows; $row_counter ++ ){
printf( "%.1f - %.1f : ", ($row_counter/$hist_rows), ((1/$hist_rows) + ($row_counter/$hist_rows)));
for( my $counter = 0; $counter < ($histogram[ $row_counter ]*(30/($n/10))); $counter ++ ){
# the extra math on the condition is to normalize the rows of the
# histogram to be 30 chars long
print "*";
}
print "\n";
}</lang>
Usage:
perl rand_statistics.pl (number of values)
$ perl rand_statistics.pl 100 100 numbers Mean: 0.531591369804339 Stddev: 0.28440375340793 0.0 - 0.1 : *************************** 0.1 - 0.2 : ************************ 0.2 - 0.3 : *************************** 0.3 - 0.4 : ************************ 0.4 - 0.5 : ********************************* 0.5 - 0.6 : ************************************ 0.6 - 0.7 : ************************************ 0.7 - 0.8 : ****************** 0.8 - 0.9 : *************************************** 0.9 - 1.0 : ************************************ $ perl rand_statistics.pl 1000 1000 numbers Mean: 0.51011452684812 Stddev: 0.29490201218115 0.0 - 0.1 : ****************************** 0.1 - 0.2 : ******************************* 0.2 - 0.3 : *************************** 0.3 - 0.4 : ***************************** 0.4 - 0.5 : ********************************** 0.5 - 0.6 : **************************** 0.6 - 0.7 : ************************ 0.7 - 0.8 : ************************************* 0.8 - 0.9 : ******************************** 0.9 - 1.0 : ********************************* $ perl rand_statistics.pl 10000 10000 numbers Mean: 0.495329167703333 Stddev: 0.285944419431566 0.0 - 0.1 : ***************************** 0.1 - 0.2 : ******************************* 0.2 - 0.3 : ********************************* 0.3 - 0.4 : ******************************* 0.4 - 0.5 : ****************************** 0.5 - 0.6 : ******************************* 0.6 - 0.7 : ****************************** 0.7 - 0.8 : ****************************** 0.8 - 0.9 : ***************************** 0.9 - 1.0 : ****************************** $ perl rand_statistics.pl 10000000 10000000 numbers Mean: 0.499973935749229 Stddev: 0.2887231680817 0.0 - 0.1 : ****************************** 0.1 - 0.2 : ******************************* 0.2 - 0.3 : ****************************** 0.3 - 0.4 : ******************************* 0.4 - 0.5 : ****************************** 0.5 - 0.6 : ******************************* 0.6 - 0.7 : ****************************** 0.7 - 0.8 : ****************************** 0.8 - 0.9 : ******************************* 0.9 - 1.0 : *******************************
PicoLisp
The following has no limit on the number of samples. The 'statistics' function accepts an executable body 'Prg', which it calls repeatedly to get the samples. <lang PicoLisp>(scl 6)
(de statistics (Cnt . Prg)
(prinl Cnt " numbers")
(let (Sum 0 Sqr 0 Hist (need 10 NIL 0))
(do Cnt
(let N (run Prg 1) # Get next number
(inc 'Sum N)
(inc 'Sqr (*/ N N 1.0))
(inc (nth Hist (inc (/ N 0.1)))) ) )
(let M (*/ Sum Cnt)
(prinl "Mean: " (round M))
(prinl "StdDev: "
(round
(sqrt
(* 1.0
(- (*/ Sqr Cnt) (*/ M M 1.0)) ) ) ) ) )
(for (I . H) Hist
(prin (format I 1) " ")
(do (*/ H 400 Cnt) (prin '=))
(prinl) ) ) )</lang>
Test: <lang PicoLisp>(statistics 100
(rand 0 (dec 1.0)) )
(prinl)
(statistics 10000
(rand 0 (dec 1.0)) )
(prinl)
(statistics 1000000
(rand 0 (dec 1.0)) )
(prinl)</lang> Output:
100 numbers Mean: 0.501 StdDev: 0.284 0.1 ======================================== 0.2 ==================================== 0.3 ==================================================== 0.4 ======================== 0.5 ======================== 0.6 ================================================================ 0.7 ======================================================== 0.8 ==================================== 0.9 ======================== 1.0 ============================================ 10000 numbers Mean: 0.501 StdDev: 0.288 0.1 ======================================= 0.2 ======================================== 0.3 ======================================= 0.4 ========================================= 0.5 ========================================= 0.6 ======================================== 0.7 ========================================= 0.8 ======================================== 0.9 ======================================== 1.0 ======================================== 1000000 numbers Mean: 0.500 StdDev: 0.289 0.1 ======================================== 0.2 ======================================== 0.3 ======================================== 0.4 ======================================== 0.5 ======================================== 0.6 ======================================== 0.7 ======================================== 0.8 ======================================== 0.9 ======================================== 1.0 ========================================
Python
The second function, sd2 only needs to go once through the numbers and so can more efficiently handle large streams of numbers. <lang python>def sd1(numbers):
if numbers:
mean = sum(numbers) / len(numbers)
sd = (sum((n - mean)**2 for n in numbers) / len(numbers))**0.5
return sd, mean
else:
return 0, 0
def sd2(numbers):
if numbers:
sx = sxx = n = 0
for x in numbers:
sx += x
sxx += x*x
n += 1
sd = (n * sxx - sx*sx)**0.5 / n
return sd, sx / n
else:
return 0, 0
def histogram(numbers):
h = [0] * 10
maxwidth = 50 # characters
for n in numbers:
h[int(n*10)] += 1
mx = max(h)
print()
for n, i in enumerate(h):
print('%3.1f: %s' % (n / 10, '+' * int(i / mx * maxwidth)))
print()
if __name__ == '__main__':
import random
for i in range(1,8):
n = [random.random() for i in range(10**i)]
print("\n##\n## %i numbers\n##" % 10**i)
print(' Naive method: sd: %8.6f, mean: %8.6f' % sd1(n))
print(' Second method: sd: %8.6f, mean: %8.6f' % sd2(n))
histogram(n)</lang>
- Section of output
for larger sets of random numbers, the distribution of numbers between the bins of the histogram evens out.
... ## ## 100 numbers ## Naive method: sd: 0.288911, mean: 0.508686 Second method: sd: 0.288911, mean: 0.508686 0.0: +++++++++++++++++++++++++++++++ 0.1: ++++++++++++++++++++++++++++ 0.2: +++++++++++++++++++++++++ 0.3: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.4: ++++++++++++++++++ 0.5: +++++++++++++++++++++++++++++++ 0.6: ++++++++++++++++++ 0.7: +++++++++++++++++++++++++++++++++++++ 0.8: ++++++++++++++++++++++++++++++++++++++++ 0.9: +++++++++++++++++++++++++++++++ ... ## ## 10000000 numbers ## Naive method: sd: 0.288750, mean: 0.499839 Second method: sd: 0.288750, mean: 0.499839 0.0: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.1: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.2: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.3: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.4: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.5: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.6: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.7: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.8: +++++++++++++++++++++++++++++++++++++++++++++++++ 0.9: +++++++++++++++++++++++++++++++++++++++++++++++++