## Statistics

### September 27, 2011

In today’s exercise we calculate some of the basic measures in statistics: mean, standard deviation, linear regression, and correlation. The only hard part is that different sources use different standard names to refer to the different statistics. The formulas are shown below; all the summations are over $i$ from 1 to the number of items $n$:

mean: $\mu = \bar{x} = \frac{1}{n} \sum x_i$

standard deviation: $\sigma = s = \sqrt{\frac{1}{n} \sum (x_i - \mu)^2}$

linear regression: $y = mx+b = \hat{\beta}x + \hat{\alpha}$

slope: $m = \hat{\beta} = \frac{n \sum x_i y_i - \sum x_i \sum y_i}{n \sum x_i^2 - \left(\sum x_i\right)^2}$

intercept: $b = \hat{\alpha} = \frac{1}{n} \sum y - \hat{\beta} \frac{1}{n} \sum x$

correlation: $r = \frac{\sum (x_i - \bar{x}) (y_i - \bar{y})}{(n-1) s_x s_y}$

Your task is to write functions to compute these basic statistics. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 8 Responses to “Statistics”

1. DGel said

The implementation of your standard deviation (and thus correlation) is wrong, given the definitions on page 1. Your definition says to divide by N, you divide by N – 1…

My implementation in Go:

 

 package main import ( "fmt" "math" ) func mean(data []float64) float64 { var sum float64 for _, x := range data {sum += x} return sum / float64(len(data)) } func sd(data []float64) float64 { mn := mean(data) var sd float64 for _, x := range data {diff := x - mn; sd += diff * diff} return math.Sqrt(sd / float64(len(data))) } func regress(xs, ys []float64) (slope, intercept float64) { var sum_xs, sum_xs2, sum_ys, sum_xs_ys, n float64 for i, _ := range xs { sum_xs += xs[i] sum_xs2 += xs[i] * xs[i] sum_ys += ys[i] sum_xs_ys += xs[i]*ys[i] } n = float64(len(xs)) slope = (n*sum_xs_ys - sum_xs*sum_ys) / (n*sum_xs2 - sum_xs*sum_xs) intercept = sum_ys/n - slope*sum_xs/n return } func correlation(xs, ys []float64) float64 { x_bar := mean(xs) y_bar := mean(ys) var cor float64 for i, _ := range xs {cor += (xs[i] - x_bar) * (ys[i] - y_bar)} return cor / (float64(len(xs) - 1)*sd(xs)*sd(ys)) } func main() { xs := []float64{1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83} ys := []float64{52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46} fmt.Println(mean(xs), mean(ys)) fmt.Println(sd(xs), sd(ys)) fmt.Println(regress(xs, ys)) fmt.Println(correlation(xs, ys)) return } 

2. Jussi Piitulainen said

I think I was taught to divide by n – 1 when the deviation from the (unknown) population mean is wanted but the (known) sample mean is used instead in the formula. The sample values are said to lose one “degree of freedom” because they can not all deviate freely from their own mean.

3. Paul Hofstra said

See: http://en.wikipedia.org/wiki/Standard_deviation
If you divide by n, the standard deviation is biased. Dividing by n-1 gives an unbiased standard deviation.

4. Jussi Piitulainen said

By way of conversation, here is an approach I find much fun. I lift constants to be vecs (indexed sequences) so that everything is uniform, and then I map binary or unary operations on these vecs. Like in the language of R but more rigidly and in Scheme. The goal is a special language that allows to explore descriptions like “the mean square deviation from the mean” in the code itself. Someone should write The Structure and Interpretation of Statistical, er, Something.

Ok, I get carried away. A variation on the theme anyway. I’ve included one of Anscombe’s cases.

(define (vec . args) (case-lambda ((k) (list-ref args k))
(() (length args))))

(define fun (case-lambda
((f u w) (case-lambda ((k) (f (u k) (w k))) (() (u))))
((f v) (case-lambda ((k) (f (v k))) (() (v))))))
I
(define (con v c) (case-lambda ((k) c) (() (v))))

(define (sum v) (con v (do ((k 0 (+ k 1)) (s 0 (+ s (v k))))
((= k (v)) s))))
(define (mean v) (fun / (sum v) (con v (v))))
(define (mean1 v) (fun / (sum v) (con v (- (v) 1))))

(define (var1 v) (let* ((d (fun - v (mean v)))) (mean1 (fun * d d))))
(define (stddev1 v) (fun sqrt (var1 v)))

(define (std1 v) (fun / (fun - v (mean v)) (stddev1 v)))
(define (cor1 u w) (mean1 (fun * (std1 u) (std1 w))))

(define xs-1 (vec 10 8 13 9 11 14 6 4 12 7 5))
(define ys-1 (vec 8.04 6.95 7.58 8.81 8.33 9.96 7.24 4.26 10.84 4.82 5.68))

(define (test)
((means ,((mean xs-1) 0) ,((mean ys-1) 0))
(devs ,((stddev1 xs-1) 0) ,((stddev1 ys-1) 0))
(corr ,((cor1 xs-1 ys-1) 0))))

;((means 9 7.500909090909093)
;  (devs 3.3166247903554 2.031568135925815)
;  (corr 0.8164205163448399))

5. DGel said

@Jussi Piitulainen, Paul Hofstra:

Yes, but that’s not how he defined standard deviation on page 1. Thus the confusion..

6. Jussi Piitulainen said

@DGel: Yes. It may be better to deviate from the definition on page 1, especially when even the model implementation does so.

7. Axio said

How ugly :)

let sum = List.fold_left (+.) 0.;;

let mean =
function
| [] -> 0.
| l -> sum l /. (float_of_int (List.length l));;

let std =
function
| [] -> 0.
| l ->
let mu = mean l in
let le = float_of_int (List.length l) in
sqrt
((List.fold_left (fun r e -> r +. ((e -. mu) *. (e -. mu))) 0. l)
/. le);;

let slope xs ys =
let le = float_of_int (List.length xs) in
let sumxs = sum xs in
let sumys = sum ys in
((le *.
List.fold_left2 (fun r xi yi -> r +. (xi *. yi)) 0. xs ys)
-.
(sumxs *.  sumys))
/.
((le *.
(List.fold_left (fun r xi -> r +. (xi *. xi)) 0. xs))
-. (sumxs *. sumxs));;

let intercept xs ys =
let le = float_of_int (List.length xs) in
((sum ys) /. le) -.
(((slope xs ys) *. (sum xs)) /. le);;

let lin_reg xs ys = ((slope xs ys), (intercept xs ys));;

let correlation xs ys =
let le = float_of_int (List.length xs) in
let mx = mean xs in
let my = mean ys in
let sx = std xs in
let sy = std ys in
(List.fold_left2 (fun r xi yi -> r +. ((xi -. mx) *.(yi -. my))) 0. xs ys)
/.
((le -. 1.) *. sx *. sy)
;;

let xs = [1.47;1.50;1.52;1.55;1.57;1.60;1.63;1.65;1.68;1.70;1.73;1.75;1.78;1.80;1.83];;
let ys = [52.21;53.12;54.48;55.84;57.20;58.57;59.93;61.29;63.11;64.47;66.28;68.10;69.92;72.19;74.46];;
let wrap point v = let eps = 0.001 in point -. eps < v && v < point +. eps;;
assert (wrap 1.6506666666666665 (mean xs));;
assert (wrap 62.078 (mean ys));;
assert (wrap 0.110 (std xs));;
assert (wrap 6.798 (std ys));;
assert (wrap 61.272186542107434 (slope xs ys));;
assert (wrap (-39.06195591883866) (intercept xs ys));;
assert (wrap 1.06562549311809596 (correlation xs ys));;
`