Triple birthday problem.ipynb

triple-birthday-problem.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Triple birthday problem.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/LiberiFatali/9660278f6a513018eba3b2da03636e71/triple-birthday-problem.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a354Fzf-CGOv",
        "colab_type": "code",
        "outputId": "5ad04dbe-75ee-4a53-a5b7-2d43241cfa5e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 33
        }
      },
      "source": [
        "### Implement the last post of http://mathforum.org/library/drmath/view/56650.html\n",
        "\n",
        "from scipy.special import comb, perm\n",
        "\n",
        "\n",
        "def prob_trio(N):\n",
        "  n_element = int(1 + N/2)      # [0 --> N/2]\n",
        "  list_pn = [0] * n_element\n",
        "  sum_list_pn = 0\n",
        "  for K in range(n_element):\n",
        "      list_pn[K] = comb(365, K) * perm(N, 2*K) * perm(365-K, N-2*K) / (2**K * 365**N)\n",
        "      sum_list_pn += list_pn[K]\n",
        "  # print(list_pn)\n",
        "  return 1-sum_list_pn\n",
        "\n",
        "N = 70\n",
        "print('Prob of at least trio: ', prob_trio(N))\n"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Prob of at least trio:  0.3064873269820343\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zVzCC8p4Cgtn",
        "colab_type": "code",
        "outputId": "cd5ddf83-d8ae-4991-c4b3-0c53992706f5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 33
        }
      },
      "source": [
        "### Take leap years into consideration\n",
        "### Implement https://drive.google.com/open?id=1c3ATed3-GhSDMSJ83S4l3ey_SNDsBNzp\n",
        "\n",
        "from scipy.special import comb, perm\n",
        "\n",
        "def prob_trio(N):\n",
        "  n_element = int(1 + N/2)      # [0 --> N/2]\n",
        "  list_pn = [0] * n_element\n",
        "  sum_list_pn = 0\n",
        "  for K in range(n_element):\n",
        "      list_pn[K] = comb(365, K) * perm(N, 2*K) * perm(365-K, N-2*K) / (2**K * 365**N)\n",
        "      sum_list_pn += list_pn[K]\n",
        "  # print(list_pn)\n",
        "  return 1-sum_list_pn\n",
        "\n",
        "def prob_feb29_exact(m_people, N):\n",
        "  p1 = 8 / (365*31 + 8)      # [1970 --> 2000]\n",
        "  p2 = comb(N, m_people)*p1**m_people * (1-p1)**(N-m_people)\n",
        "  return p2\n",
        "\n",
        "def prob_feb29_atleast(m_people, N):\n",
        "  p1 = 8 / (365*31 + 8)      # [1970 --> 2000]\n",
        "  p2 = comb(N, m_people)*p1**m_people\n",
        "  return p2\n",
        "\n",
        "N = 70\n",
        "final_prob = prob_feb29_exact(0, N)*prob_trio(N) + prob_feb29_exact(1, N)*prob_trio(N-1) + prob_feb29_exact(2, N)*prob_trio(N-2) + prob_feb29_atleast(3, N)\n",
        "print('Include leap years. Prob of at least trio: ', final_prob)\n"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Include leap years. Prob of at least trio:  0.30597723013004147\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}