The Iron Bar
October 6, 2015
Our experiment will perform n trials, each of a list of 1001 numbers from 1 to 100, comparing the iron-bar median to the actual median determined by sorting:
(define (rand-list n)
(let loop ((n n) (xs (list)))
(if (zero? n) xs
(loop (- n 1) (cons (+ (randint 100) 1) xs)))))
(define (iron-bar n)
(do ((i 1 (+ i 1))) ((< n i))
(let* ((xs (rand-list 1001)) (actual-median (list-ref (sort < xs) 500)))
(let loop ((med (car xs)) (xs (cdr xs)))
(if (null? xs)
(for-each display `(,i ": actual median = "
,actual-median ", iron-bar median = " ,med #\newline))
(loop (+ (if (< (car xs) med) -1 1) med) (cdr xs)))))))
And here are some tests of the iron bar:
> (iron-bar 10) 1: actual median = 48, iron-bar median = 49 2: actual median = 53, iron-bar median = 51 3: actual median = 51, iron-bar median = 50 4: actual median = 51, iron-bar median = 50 5: actual median = 52, iron-bar median = 52 6: actual median = 54, iron-bar median = 53 7: actual median = 51, iron-bar median = 54 8: actual median = 50, iron-bar median = 59 9: actual median = 52, iron-bar median = 48 10: actual median = 51, iron-bar median = 47
That’s not bad. Most of the results are within 1 or 2 of the correct answer. Trial 8 is quite bad: that happens when there is a long streak near the end of the data.
We used the random number generator in the Standard Prelude. You can run the program at http://ideone.com/jD90tT.
Haskell:
Scheming with Gambit and a biased distribution. Stealing signum from from above.
Thinking the result seems more or less as it probably should.
fun median [] = ~1 | median (x::xs) = let fun loop med [] = med | loop med (x::xs) = if x = med then loop med xs else if x > med then loop (med+1) xs else loop (med-1) xs in loop x xs end val rand_list = [ 3,5,68,69,34,48,56,2,65,28,20,77,13,27,58,0,70,61,36,54,85,1,29,42,87,12,49,58,1,15,62,29,74,32,91,25,9,48,94,2,21,7,66,18,11,22,3,72,87,48,65,71,15,62,85,13,28,55,9,22,88,91,65,85,30,40,56,92,77,41,76,23,98,69,55,7,89,4,27,33,44,31,2,36,98,84,32,80,22,96,55,22,45,10,54,38,55,1,40,17,44,81,17,49,57,28,33,15,11,44,78,69,56,74,1,4,31,91,13,58,55,86,68,82,82,40,15,93,38,36,11,89,11,8,64,85,69,25,87,57,46,26,75,4,50,77,90,59,76,29,82,16,61,53,30,34,98,60,58,79,70,30,82,38,85,85,42,42,68,46,70,48,44,41,56,42,12,31,38,69,90,89,64,93,10,13,93,38,59,0,87,35,9,37,53,61,97,20,76,96,40,41,76,98,84,76,44,11,37,78,25,4,1,46,86,26,67,13,16,62,77,2,27,21,80,16,76,38,83,42,14,85,8,97,0,97,44,73,94,47,39,46,71,46,51,76,60,24,53,40,18,87,32,75,51,97,50,65,5,65,64,60,16,91,44,8,52,44,39,83,60,84,52,85,66,69,81,83,41,70,88,67,5,72,90,71,49,28,26,80,73,68,97,78,55,57,86,96,50,65,79,82,35,84,92,74,59,3,99,98,31,67,29,63,3,34,68,51,4,57,35,20,49,51,20,91,78,88,83,29,65,70,4,93,15,70,12,36,46,67,81,1,58,7,89,36,72,0,92,79,10,3,30,87,29,18,66,40,74,89,33,1,59,4,32,69,44,2,23,20,23,99,48,7,5,84,63,47,27,71,5,94,93,35,67,15,3,23,80,48,11,81,62,97,57,41,82,71,71,51,96,93,32,38,92,76,9,37,13,53,29,50,66,21,80,58,19,86,53,0,97,9,83,62,90,19,19,63,15,16,86,24,41,92,74,14,65,85,13,32,82,26,0,89,95,61,65,83,59,33,33,22,26,63,9,45,38,7,56,84,51,20,87,81,64,1,17,89,20,26,92,70,60,89,97,4,20,5,26,10,74,84,76,36,93,27,58,32,23,17,29,64,1,9,8,1,93,73,64,1,20,72,22,44,67,37,9,1,64,45,40,88,53,99,84,33,91,6,23,76,96,88,3,66,71,27,81,50,98,79,28,63,95,58,33,15,3,81,18,73,27,1,47,0,7,34,44,86,93,59,13,74,65,74,82,1,90,7,65,87,67,18,57,64,65,37,97,36,14,65,82,8,15,42,11,9,56,97,33,26,64,55,3,48,59,62,97,14,2,43,4,75,77,30,68,26,15,48,4,0,16,7,80,66,69,26,34,42,58,52,38,91,44,71,0,59,82,52,76,67,88,25,56,37,26,84,34,61,35,62,78,36,41,25,80,72,35,38,66,31,98,92,75,8,25,91,62,5,8,30,67,51,51,74,48,11,61,75,78,85,7,39,69,16,39,71,69,77,26,42,61,90,21,67,19,47,96,40,38,54,61,92,47,3,25,41,99,43,61,76,18,33,70,59,95,35,3,62,8,27,42,30,40,51,25,89,28,12,98,18,64,65,39,75,72,55,51,5,76,94,18,18,23,29,16,70,56,80,55,18,49,11,71,91,92,34,95,81,94,2,22,82,68,40,25,31,28,74,83,35,16,76,75,20,20,58,74,5,96,74,43,83,66,10,51,68,85,83,7,12,65,7,38,8,42,62,39,3,18,49,52,51,77,2,63,35,7,93,21,21,48,27,6,74,67,26,8,44,18,23,57,68,14,77,61,15,23,68,78,55,39,66,63,49,58,68,88,64,22,27,84,5,10,42,58,1,6,15,65,15,73,94,75,13,71,28,39,69,33,4,93,97,11,17,50,26,46,46,27,78,13,51,70,97,13,34,23,13,79,98,21,84,20,6,60,88,94,14,92,37,35,24,64,66,90,16,67,69,39,13,67,4,70,47,76,86,25,10,1,22,27,77,94,48,16,60,89,98,21,41,62,40,46,59,5,1,28,36,33,43,88,84,51,21,59,27,45,22,40,55,86,47,7,21,53,7,56,5,62,34,50,47,23,35,64,36,62,77,21,57,10,41,14,91,80,29,6,75,88,3,1,83,2,40,90,76,94,65,19,49,16,86,95,23,95,12,52,67,55,30,11,40,51,36,52,59,6,16,94,44,54,69,16,91,7,35,14,22,15,52,70,11,72,95,79,99,19,16,43,90,90 ]; val estimated_median = median rand_list; (* 47 *) (* actual median 49 *)Maybe I need more coffee and to read more carefully, but isn’t this just a Beta-Bernoulli setup from Bayesian statistics?
@Graham: I rather thought the whole article was a joke. I actually checked to be sure it wasn’t published April 1 and somebody reposted it late. They even use jargon: volume-conserving, TOW dynamics, quantum dots. And they really lost me when they claimed their iron bar can compute the probabilities of winning on a slot machine that it doesn’t even play. At least I managed to get a decent exercise out of it, and maybe even amuse some of my readers.
A one-liner in J:
> median =: [ + (* @: -)~
Another Haskell version, somewhat golfed.
#include <vector> #include <numeric> #include <algorithm> #include <iostream> #include <random> #include <iterator> void approximate_median(const std::vector<int>& numbers) { if (numbers.empty()) return; double real_median = std::accumulate(numbers.begin(), numbers.end(), 0.0) / numbers.size(); int approximation = std::accumulate( numbers.begin(), numbers.end(), numbers.front(), [](auto median, auto n){ if (n > median) return ++median; else if (n < median) return --median; else return median; }); std::cout << "Approximate median = " << approximation << "; actual median = " << real_median << '\n'; } int main() { using E = std::default_random_engine; using D = std::uniform_int_distribution<int>; auto generator = [e = E(), d = D(1, 100)]() mutable { return d(e); }; std::vector<int> numbers; std::generate_n(std::back_inserter(numbers), 1000, std::ref(generator)); approximate_median(numbers); }Dang. Confused median with mean. Second try:
#include <vector> #include <numeric> #include <algorithm> #include <iostream> #include <random> #include <iterator> int median(const std::vector<int>& numbers) { std::vector<int> ns(numbers.begin(), numbers.end()); auto halfway = ns.size() / 2; std::nth_element(ns.begin(), ns.begin() + halfway, ns.end()); if (ns.size() % 2) { return ns[halfway]; } else { return (ns[halfway - 1] + ns[halfway]) / 2; } } int approximate_median(const std::vector<int>& numbers) { if (numbers.empty()) return 0; return std::accumulate(numbers.begin(), numbers.end(), numbers.front(), [](auto median, auto n){ if (n > median) return ++median; else if (n < median) return --median; else return median; }); } int main() { using E = std::default_random_engine; using D = std::uniform_int_distribution<int>; auto generator = [e = E(), d = D(1, 100)]() mutable { return d(e); }; std::vector<int> numbers; std::generate_n(std::back_inserter(numbers), 1000, std::ref(generator)); auto real_median = median(numbers); auto approximation = approximate_median(numbers); std::cout << "Real median = " << real_median << "; approximation = " << approximation << '\n'; }