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.

Pages: 1 2

9 Responses to “The Iron Bar”

  1. Francesco said

    Haskell:

    m []     jn       = maybe 0 id jn
    m (a:as) Nothing  = m as (Just a)
    m (a:as) (Just k) = m as (Just $ k + (signum $ a-k))
    
  2. Jussi Piitulainen said

    Scheming with Gambit and a biased distribution. Stealing signum from from above.

    (define (sign x) (if (negative? x) -1 (if (positive? x) 1 0)))
    (define (fold f x xs) (if (pair? xs) (fold f (f x (car xs)) (cdr xs)) x))
    (define (step x y) (+ x (sign (- y x))))
    (define (walk list) (fold step (car list) (cdr list)))
    (define (skew) (abs (- (random-integer 100) (random-integer 50))))
    (define (make n) (do ((k 0 (+ k 1)) (data '() (cons (skew) data))) ((= k n) data)))
    ;;; Displays a list of approximate medians for 10 samples of 100
    ;;; numbers between 0 (inclusive) and 100 (exclusive) with a heavy
    ;;; bias towards the smaller numbers.
    (display (do ((k 0 (+ k 1)) (meds '() (cons (walk (make 100)) meds))) ((= k 10) meds)))
    (newline)
    

    Thinking the result seems more or less as it probably should.

    $ gsi med.scm
    (29 26 31 30 30 28 28 31 24 32)
    
  3. mcmillhj said
    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 *)
    
    
  4. Graham said

    Maybe I need more coffee and to read more carefully, but isn’t this just a Beta-Bernoulli setup from Bayesian statistics?

  5. programmingpraxis said

    @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.

  6. agambrahma said

    A one-liner in J:

    > median =: [ + (* @: -)~

  7. Globules said

    Another Haskell version, somewhat golfed.

    import Data.List (foldl1')
    
    -- A running estimate of the median of a list.  (Golfed version.)
    med :: [Int] -> Int
    med [] = 0
    med xs = foldl1' (\m x -> fromEnum (x `compare` m) + m - 1) xs
    
  8. fisherro said
    #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);
    }
    
  9. fisherro said

    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';
    }
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: