Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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"
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment