[c++]Monty-Hall simulation

mh.cpp

// clang++ -Wall -Wextra -pedantic -std=c++17 mh.cpp -o mh
#include <vector>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <random>

std::default_random_engine engine;

// utilities
template<typename I0, typename I1> I0 randint(I1 start, I0 end) {
  return std::uniform_int_distribution<I0>(start, end - 1)(::engine);
}
template<typename Vec, typename E> bool contains(Vec& vec, E& elem) {
  return std::find(vec.begin(), vec.end(), elem) != vec.end();
}

// fake c++20 range/view
template<typename Vec, typename View> Vec operator|(Vec a, View v) {
  return v(a);
};
template<typename I> auto iota(I n) {
  return [n](auto a) {
    std::iota(a.begin(), a.end(), n);
    return a;
  };
}
template<typename Pred> auto filter(Pred p) {
  return [p](auto a) { 
    decltype(a) r(0);
    std::copy_if(a.begin(), a.end(), std::back_inserter(r), p);
    return r;
  };
}
template<typename I> auto take(I n) {
  return [n](auto a) {
    decltype(a) r(n);
    std::copy(a.begin(), a.begin() + n, r.begin());
    return r;
  };
}
template<typename MapF> auto transform(MapF m) {
  return [m](auto a) { 
    decltype(a) r(a.size());
    std::transform(a.begin(), a.end(), r.begin(), m);
    return r;
  };
}

// randomly choosen n-elements from vector a
template<typename Vec, typename I> Vec choose(Vec a, I n) {
  assert(a.size() >= decltype(a.size())(n));
  if (a.size() == 0) return Vec(0);
  auto idx = Vec(a.size()) | ::iota(0);
  for (auto i = 0; i < n; ++i) {
    auto j = ::randint(i, a.size());
    std::iter_swap(idx.begin() + i, idx.begin() + j);
  }
  return idx | ::take(n) | ::transform([&](auto i){return a[i];});
}
template<typename I> auto choose(I n) {
  return [n](auto a){return choose(a, n);};
}


// no re-choice
bool mh0(int doors, int cars, int opens) {
  assert(doors > cars + opens);
  auto box = std::vector<int>(doors) | ::iota(0);
  auto carIndex = box | ::choose(cars);

  auto choice1 = ::randint(0, doors);
  return ::contains(carIndex, choice1);
}

// MC opens random doors then re-choice randomly
bool mh1(int doors, int cars, int opens) {
  assert(doors > cars + opens);
  auto box = std::vector<int>(doors) | ::iota(0);
  auto carIndex = box | ::choose(cars);

  auto choice1 = ::randint(0, doors);

  auto mcBox = box | ::filter([&](auto i){
    return i != choice1 && !::contains(carIndex, i);
  });
  auto openIndex = mcBox | ::choose(opens);

  auto chanceBox = box | ::filter([&](auto i){
    return i != choice1 && !::contains(openIndex, i);
  });
  auto choice2 = (chanceBox | ::choose(1))[0];

  return ::contains(carIndex, choice2);
}

// MC opens from earlier doors then re-choice the earliest door
// from https://qiita.com/legohasiri/items/15741d67895f4e0a594c
bool mh2(int doors, int cars, int opens) {
  assert(doors > cars + opens);
  auto box = std::vector<int>(doors) | ::iota(0);
  auto carIndex = box | ::choose(cars);

  auto choice1 = ::randint(0, doors);

  auto mcBox = box | ::filter([&](auto i){
    return i != choice1 && !::contains(carIndex, i);
  });
  std::vector<int> openIndex(mcBox.begin(), mcBox.begin() + opens);

  auto chanceBox = box | ::filter([&](auto i){
    return i != choice1 && !::contains(openIndex, i);
  });
  auto choice2 = chanceBox[0];

  return ::contains(carIndex, choice2);
}


// simulation and output
void sim(int n, int doors, int cars, int opens) {
  std::cout << "[doors " << doors << ", cars " << cars <<
    ", opens " << opens << "]" << std::endl;
  int s0 = 0, s1 = 0, s2 = 0;
  for (auto i = 0; i < n; ++i) s0 += ::mh0(doors, cars, opens);
  for (auto i = 0; i < n; ++i) s1 += ::mh1(doors, cars, opens);
  for (auto i = 0; i < n; ++i) s2 += ::mh2(doors, cars, opens);
  std::cout << "No rechoice: " << double(s0) / n * 100 << "%" << std::endl;
  std::cout << "MC random  : " << double(s1) / n * 100 << "%" << std::endl;
  std::cout << "MC earlier : " << double(s2) / n * 100 << "%" << std::endl;
  std::cout << std::endl;
}


// output
int main(void) {
  int n = 10000;

  {
    int doors = 3, cars = 1, opens = 1;
    sim(n, doors, cars, opens);
    // No rechoice: 1/3
    // MC random  : 1/3*0/1 + 2/3*1/1 = 2/3
    // MC eralier : 1/3*0/1 + 2/3*1/1 = 2/3
  }
  {
    int doors = 4, cars = 1, opens = 1;
    sim(n, doors, cars, opens);
    // No rechoice: 1/4
    // MC random  : 1/4*0/2 + 3/4*1/2 = 3/8 = 37.5%
    // MC eralier : 1/4*0/2 + 3/4*(1/3 + 2/3*1/2) = 1/2
  }
  {
    int doors = 4, cars = 2, opens = 1;
    sim(n, doors, cars, opens);
    // No rechoice: 2/4
    // MC random  : 2/4*1/2 + 2/4*2/2 = 3/4
    // MC eralier : 2/4*(1/3 + 2/3*1/2) + 2/4*2/2 = 5/6 = 83.33%
  }
  {
    int doors = 5, cars = 1, opens = 1;
    sim(n, doors, cars, opens);
    // No rechoice: 1/5
    // MC random  : 1/5*0/3 + 4/5*1/3 = 4/15 = 26.66%
    // MC eralier : 1/5*0/3 + 4/5*(1/4 + 3/4*1/3) = 2/5
  }
  {
    int doors = 5, cars = 2, opens = 1;
    sim(n, doors, cars, opens);
    // No rechoice: 2/5
    // MC random  : 2/5*1/3 + 3/5*2/3 = 8/15 = 53.33%
    // MC eralier : 2/5*(1/4+3/4*1/3) + 3/5*(2/4 + 2/4*(1/3+2/3*1/2)) = 7/10
  }
  {
    int doors = 5, cars = 3, opens = 1;
    sim(n, doors, cars, opens);
    // No rechoice: 3/5
    // MC random  : 3/5*2/3 + 2/5*3/3 = 4/5
    // MC eralier : 3/5*(2/4+2/4*2/3) + 2/5*(3/4 + 1/4*3/3) = 9/10
  }
}

result.txt

[doors 3, cars 1, opens 1]
No recoice: 33.51%
MC random : 67.46%
MC earlier: 66.91%

[doors 4, cars 1, opens 1]
No recoice: 25.37%
MC random : 37.41%
MC earlier: 50.34%

[doors 4, cars 2, opens 1]
No recoice: 49.48%
MC random : 74.56%
MC earlier: 83.63%

[doors 5, cars 1, opens 1]
No recoice: 19.8%
MC random : 26.69%
MC earlier: 40.59%

[doors 5, cars 2, opens 1]
No recoice: 39.88%
MC random : 53.4%
MC earlier: 70.51%

[doors 5, cars 3, opens 1]
No recoice: 59.75%
MC random : 80.01%
MC earlier: 89.93%

[Probabilities]

MC random 4-door 1-car case: 1/4 * 0/2 + 3/4 * 1/2 = 3/8
MC earlier 4-door 1-car case: 1/4 * 0/2 + 3/4 * (1/3 + 2/3 * 1/2) = 1/2

MC random 5-door 1-car case: 1/5 * 0/3 + 4/5 * 1/3 = 4/15
MC earlier 5-door 1-car case: 1/5 * 0/3 + 4/5 * (1/4 + 3/4 * 1/3) = 2/5

MC random 4-door 2-car case: 2/4 * 1/2 + 2/4 * 2/2 = 3/4
MC earlier 4-door 2-car case: 2/4*(1/3+2/3*1/2) + 2/4 * 2/2 = 5/6

MC random 5-door 2-car case: 2/5 * 1/3 + 3/5 * 2/3 = 8/15
MC earlier 5-door 2-car case: 2/5*(1/4+3/4*1/3) + 3/5 * (2/4 + 2/4 * (1/3 + 2/3 * 1/2) = 1/5 + 3/10 * (1 +2/3) = 1/5 + 1/2 = 7/10

MC random 5-door 3-car case: 3/5 * 2/3 + 2/5 * 3/3 = 4/5
MC earlier 5-door 3-car case: 3/5 * (2/4 + 2/4 * 2/3) + 2/5 * (3/4 + 1/4 * 3/3) = 3/10 * (1+2/3) + 2/5 = 1/2+2/5 = 9/10