{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3)\n",
    "## e)\n",
    "The event B is defined by \"some two people in the group share a birthday\". The complement of B is: \"no two people in teh group share a birthday\". All the sequences of n birthdays that satisfy the complement of B are all the permutations  without replacement of n from 365. \n",
    "\n",
    "Permutations without replacement are given by: $_nP_k = \\frac{n!}{(n-k)!}$\n",
    "\n",
    "Then, $P(B) = 1-\\frac{_{365}P_n}{365^n}$\n",
    "\n",
    "We are looking for the smallest $n$ such that $P(B) > 0.9$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "perm_wo_rep <- function(n, k) {\n",
    "    return(choose(n=n, k=k)*factorial(k))\n",
    "}\n",
    "\n",
    "prob_b <- function(n) {\n",
    "    return (1-perm_wo_rep(365, n)/365^n)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "NaN"
      ],
      "text/latex": [
       "NaN"
      ],
      "text/markdown": [
       "NaN"
      ],
      "text/plain": [
       "[1] NaN"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prob_b(100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the issue that we cannot in fact compute `perm_wo_rep` for very large numbers directly. We would need to come up with a way to compute the factorial after simplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.903151611481736"
      ],
      "text/latex": [
       "0.903151611481736"
      ],
      "text/markdown": [
       "0.903151611481736"
      ],
      "text/plain": [
       "[1] 0.9031516"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=41\n",
    "1-perm_wo_rep(365, n)/(365^n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "66430"
      ],
      "text/latex": [
       "66430"
      ],
      "text/markdown": [
       "66430"
      ],
      "text/plain": [
       "[1] 66430"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "choose(n=365, k=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.503326508108936"
      ],
      "text/latex": [
       "0.503326508108936"
      ],
      "text/markdown": [
       "0.503326508108936"
      ],
      "text/plain": [
       "[1] 0.5033265"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=23\n",
    "1-perm_wo_rep(365, n)/(365^n)- n^2/(365^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to pick a bday\n",
    "estimate_p_c <- function(n, trials=10000) {\n",
    "    c = 0\n",
    "    for (i in 1:trials) {\n",
    "        print(i)\n",
    "        if(max(table(sample(1:365, n, replace=TRUE)))>=3) {\n",
    "            c=c+1\n",
    "        }\n",
    "    }\n",
    "    return (c/trials)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = seq(1, 200)\n",
    "\n",
    "estimates = lapply(N, estimate_p_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [],
      "text/latex": [],
      "text/markdown": [],
      "text/plain": [
       "integer(0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "which(estimates>0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.0475390156062425"
      ],
      "text/latex": [
       "0.0475390156062425"
      ],
      "text/markdown": [
       "0.0475390156062425"
      ],
      "text/plain": [
       "[1] 0.04753902"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "choose(13, 2)*(choose(4, 2)^2)*choose(11, 1)*choose(4, 1)/choose(52, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "41184"
      ],
      "text/latex": [
       "41184"
      ],
      "text/markdown": [
       "41184"
      ],
      "text/plain": [
       "[1] 41184"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "choose(13, 3)*(choose(4, 2)^2)*choose(4, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}