ECE581 HW1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3713\n"
     ]
    }
   ],
   "source": [
    "result=10\n",
    "count=0\n",
    "n=10000\n",
    "for i in range(n):\n",
    "  #  print(\"i= %d\" %(i))\n",
    "    result=10\n",
    "    while ((result>0) & (result<20)):\n",
    "        out=np.floor(np.random.uniform(0,37))\n",
    "  #      print(out)\n",
    "        if out>18 :\n",
    "            result=result+1\n",
    "        else:\n",
    "            result=result-1\n",
    "  #      print(result)\n",
    "    if result==20:\n",
    "        count=count+1\n",
    "print(count/n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in ns:\n",
    "    if sda>0:\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "when simulate 100 trails, the predicted pi is 3.076923\n",
      "when simulate 1000 trails, the predicted pi is 2.952030\n",
      "when simulate 10000 trails, the predicted pi is 3.158310\n",
      "when simulate 100000 trails, the predicted pi is 3.142924\n"
     ]
    }
   ],
   "source": [
    "d=5\n",
    "l=2\n",
    "n_set=[100,1000,10000,100000]\n",
    "for n in n_set:\n",
    "    distance=np.random.uniform(0,d/2,n)\n",
    "    angle=np.random.uniform(0,np.pi/2,n)\n",
    "    cross_num=sum(l/2*np.sin(angle)-distance>0)\n",
    "    print(\"when simulate %d trails, the predicted pi is %f\" %(n,2*l/(d*cross_num/n)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_prob(trail_times,num_people):\n",
    "    count=0\n",
    "    for i in range(trail_times):\n",
    "        birthdays=np.floor(np.random.uniform(0,366,num_people))\n",
    "        if len(np.unique(birthdays))==len(birthdays):\n",
    "            count=count+1\n",
    "    prob=count/trail_times\n",
    "    return prob\n",
    "\n",
    "def P2(n):\n",
    "    return np.exp(np.log(np.arange(1,365+1)).sum()-np.log(np.arange(1,365-n+1)).sum()-n*np.log(365))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mathematical result: $ P(n) = \\frac{365！}{365^n \\times (365-n)!}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2a353a83240>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X18VOWd9/HPLyEhICRFwCcihPYGFCGgstaV3cJ9Iy1PaqUtyuJLV+26vW/RipWqa001tLW1tVhaWutW0Vq7ai20ELAtsPWhgKtgETEIUkUJWoGICUqeIL/7j5nESTKTTJJ5zvf9emHOnHOY+Xlm+Oaa61znOubuiIhIZslKdgEiIhJ7CncRkQykcBcRyUAKdxGRDKRwFxHJQAp3EZEMpHAXEclACncRkQykcBcRyUC9kvXCgwYN8qKiomS9vIhIWtqyZctBdx/c0X5JC/eioiI2b96crJcXEUlLZvZWNPupW0ZEJAMp3EVEMpDCXUQkAyWtz11E0l9DQwMVFRXU1tYmu5SMk5eXR2FhITk5OV36+wp3EemyiooK+vfvT1FREWaW7HIyhrtTWVlJRUUFw4cP79JzqFtGRLqstraWgQMHKthjzMwYOHBgt74RKdxFpFsU7PHR3eOatuHedHvASD9FRHqyDsPdzB40s/1mtj3CdjOzJWa228y2mdlZsS+zpcVrd7H84cVU3zUK7hyALz6D5Q8v5od/2klpWTk//NPOeJcgIinqy1/+MuXl5TF5rqKiIg4ePNjuPt/5zndaPD7vvPNi8trdFU3L/SFgWjvbpwMjgn+uAX7W/bIic3eG7itj+pt3kV/3dwzHqiqY/eadXL1pCpWbfsW6He+xeK0CXiTVtP5mHY9v2r/4xS8YPXp0zJ83ktbhvnHjxoS9dns6DHd3fxZ4v51dLgJ+6QHPA58ws5NjVWBrZsbsQw/S1+pbrYcCP8y9uT/lkgM/orr2KI2NjfEqQ0Q6afHaXZSWlbfoQi0tK2fx2l1dfs6PPvqImTNnMm7cOMaMGcPjjz/O5MmTm6c26devHzfffDNnn302559/Pi+88AKTJ0/mk5/8JCtXrgTgoYceYv78+c3POWvWLJ5++uk2r/X5z3+es88+mzPOOIP7778fgFtuuYWamhrGjx/PvHnzml+z6f9v4cKFjBkzhrFjx/L4448D8PTTTzN58mS++MUvctpppzFv3ry4/JKLRZ/7EGBvyOOK4Lo2zOwaM9tsZpsPHDjQ5Re0qoqI27KAy3uto6R8Fit+eW+3PjgiEhvuTnVtA8s27GkO+NKycpZt2EN1bUOXw+0Pf/gDp5xyCi+//DLbt29n2rSWnQwfffQRkydPZsuWLfTv359vfOMbrF27lhUrVlBSUtKp13rwwQfZsmULmzdvZsmSJVRWVvLd736XPn36sHXrVh599NEW+y9fvpytW7fy8ssvs27dOhYuXMi7774LwF//+lfuvfdeysvLeeONN9iwYUOX/v/bE4twD3dKN+w75e73u/sEd58weHCHk5pF5AWFHRZkNe8z/c27GLqvTCdZRZLMzCiZNZorJxaxbMMeht+6hmUb9nDlxCJKZo3u8siQsWPHsm7dOm6++Waee+45CgoKWmzPzc1tDvyxY8cyadIkcnJyGDt2LHv27OnUay1ZsoRx48Zx7rnnsnfvXl5//fV29//LX/7C3Llzyc7O5sQTT2TSpEm8+OKLAJxzzjkUFhaSlZXF+PHjO11LNGIR7hXAqSGPC4F3YvC8Ybk7ywdcxRHP7XDfvlbP7PeWaKiWSApoCvhQ3Ql2gJEjR7JlyxbGjh3LrbfeSmlpaYvtOTk5zc+flZVF7969m5ePHj0KQK9evVp04YYbW/7000+zbt06Nm3axMsvv8yZZ57Z4Rj09hqVTXUAZGdnN9cSS7EI95XA5cFRM+cCVe7+bgyeNywz4+0hs3hq+K0cyc4P/xUhVM0h1jz6o3iVIyJRauqKCRXaB98V77zzDn379uWyyy7jpptu4qWXXur0cxQVFbF161YaGxvZu3cvL7zwQpt9qqqqGDBgAH379uW1117j+eefb96Wk5NDQ0NDm7/zmc98hscff5xjx45x4MABnn32Wc4555xO19dVHU4/YGb/BUwGBplZBfBNIAfA3e8D1gAzgN3AEeDKeBXbZMHUkbgvwOxG/vqzqxj/3m/D9g1BoIvmvLeW0th4HVlZaTusXySthfaxN3XFND2GrrfgX3nlFRYuXEhWVhY5OTn87Gc/46abburUc0ycOJHhw4czduxYxowZw1lntR3NPW3aNO677z6Ki4sZNWoU5557bvO2a665huLiYs4666wW/e4XX3wxmzZtYty4cZgZd999NyeddBKvvfZap/8/u8KS1R89YcIEj9XNOhpffpyalTfR92g14T4fDiwv+iZvD5nFgqkjY/KaIgI7duzg9NNPj2rfxWt3UV3b0BzkTYGfn5ejf5cRhDu+ZrbF3Sd09HczYuKwrHGXcP/+M7nm+fM57lhVm+0GTH/zLp5ygi1+9cGLJFrgG7c3//tr6oPXv8f4yJh+igVTR9Lnwu9Tb73Dbu9r9Ux55z7uXdf+GW4RiZ/WQa5gj5+MCXd3Z9HbY7ip7uqIJ1nz69+jurZBFzeJSMbLmHA3M/Lzchj4j5dBhHHwZnD70O0sWr1DFzeJSEbLmHCHQNfM7TNPZ/mAq2kM03w3oLrs9m5fFScikuoyKtwhcHHC20NmhR01A1BQ/x73nLZTJ3JEJKNlXLgD3HD+CKpyTgy7zQxm77sbe+U3Ca5KRGLtgw8+4Kc//SkQuIp01qxZSa4odWRkuJsZG4uujThyxhpq8PWlYbeJSPoIDfd4icfUAImQkeEOMGPeV+n1+R9Hnp6gqkJ97iKJtu0JWDwG7vhE4Oe2J7r1dLfccgt/+9vfGD9+PAsXLuTDDz8MO5Xuli1bmDRpEmeffTaf+9znmmdn3Lp1K+eeey7FxcVcfPHFHDp0CIDJkyfzH//xH0yaNIlvf/vbDB8+vHmKgerqaoqKisJOOZBKMjbcIXBx0+HeJ4XdVpPdX33uIom07QlYdT1U7QU88HPV9d0K+O9+97t86lOfYuvWrXz/+98PO5VuQ0MD1113HU8++SRbtmzhqquu4rbbbgPg8ssv53vf+x7btm1j7Nix3Hnnnc3P/cEHH/DMM8/wzW9+k8mTJ7N69WoAHnvsMb7whS+Qk5PTrcMRbxkd7gD5MxfhWW3fhD7UdrvVICKdsL4UGmparmuoCayPkXBT6e7cuZPt27czdepUxo8fz7e+9S0qKiqoqqrigw8+YNKkSQBcccUVPPvss83PdckllzQvf/nLX2bZsmUALFu2jCuvjPsUWt2WEdMPtKt4TmDemcaW0xLYsXqqV99OfvGcJBUm0sNEuslOOzff6axwU+m6O2eccQabNm1q+bJVbacqCXXcccc1L0+cOJE9e/bwzDPPcOzYMcaMGROzmuMl41vu7k6fo9Vht/Wve09Xq4okSqSb7HRw85329O/fn8OHD7e7z6hRozhw4EBzuDc0NPDqq69SUFDAgAEDeO655wB45JFHmlvx4Vx++eXMnTs3LVrt0APC3cza+VANYdHqHfzwT7qZtkjcTSmBnD4t1+X0CazvooEDBzJx4kTGjBnDwoULw+6Tm5vLk08+yc0338y4ceMYP358802sH374YRYuXEhxcTFbt25t99Z78+bN49ChQ8ydO7fL9SZSRkz526FtT+CrrsdC+vscOJJdwK0189h94nSmjj6RBVNHJaYekQzRmSl/gcB5rvWlga6YgsJAsKdJ1+iTTz7J73//ex555JGEvWaPn/K3Iz72Syx/qYLJb9zD8fYhZoGpCI47VsXduQ/w9feg+pOXtZiOVETioHhO2oR5qOuuu46nnnqKNWvWJLuUqGV8uDffAea1UWzNPw6r/7DF9jzqKO33WwpmfUfBLiJh/fjHP052CZ3WI/rc8/NyuHJiEQX1+8PuU1C/X8Eu0kW6GDA+untcM77lDoHZIhsbG6l65QQ+Uf9em+3V9KN/Y6PusSrSSXl5eVRWVjJw4EA1kGLI3amsrCQvL6/Lz9Ejwh0Cs0VuGHYtU1+/k1yOtdjWx2tY8ct7mX2FbsEn0hmFhYVUVFRw4MCBZJeScfLy8igs7Pow0Z4xWibEkW8Npe/RthcvVPc+ifxbNSRSRFJbtKNlelw/RN8IFzTl17XtrhERSVc9LtzjcZWciEiq6XnhHu4qOYD6j3BNJCYiGaLnhXvxHLhgCUey81vO9V7zPg0r5rPm0R8lqzIRkZjpeeFO4IrV+uw+tB4Xk+t1nPfWUo3bFZG01yPD3cx0QZOIZLQeGe4AFuEEalXuCSxeuyvB1YiIxFaPDffG/3M7tbS8gXYtvSn58AtU1zaoa0ZE0lpU4W5m08xsp5ntNrNbwmwfamZ/NrO/mtk2M5sR+1Jjx91Z9PYYvl5/NVX0xx3c4YjnUPiJPtw+83R1zYhIWusw3M0sG1gKTAdGA3PNbHSr3b4BPOHuZwKXAj+NdaGx1DSZ2KSRg8nPORqYAtjgePuQBbVLydr+pFruIpLWomm5nwPsdvc33L0eeAy4qNU+DuQHlwuAd2JXYnwsmDqS2YceaHEDD4Ccxlp8/Z2UlpWr711E0lY04T4E2BvyuCK4LtQdwGVmVgGsAa6LSXVx5O5QtS/8tqp9LNuwR33vIpK2ogn3cJ3PrRNvLvCQuxcCM4BHzKzNc5vZNWa22cw2J3sWOTPjcO8Tw257p3EgV04sUt+7iKStaMK9Ajg15HEhbbtdrgaeAHD3TUAeMKj1E7n7/e4+wd0nDB48uGsVx1D+zEV4q6kI3KGv1XL70O0sWr1DXTMikpaiCfcXgRFmNtzMcgmcMF3Zap+3gSkAZnY6gXBP/Qmei+dgFyzB+wxo/irSdGK1fsV1VG76lbpmRCQtdRju7n4UmA/8EdhBYFTMq2ZWamYXBnf7GvBvZvYy8F/Av3qaJKKP/RJVx3Lb9D3lUcdteb9R14yIpKWo7sTk7msInCgNXVcSslwOTIxtaYnR3lQEJzQexHTrPRFJQ0ou2p+KIE2+gIiItKBwB3xKCfXWciqCegtMRVBaVq6AF5G0o3AHrHgO6/7XbXwUMsd7Tt5xgStY83LU5y4iaSeqPveeYMbYk/E9DdixwGOreZ/Z++7GzioERia1NhGRzlLLvcn60jZTEVhDDawvTVJBIiJdp3BvUlXRufUiIilM4d4kwoiZiOtFRFKYwr3JlBJoNRUBOX0C60VE0ozCvUnxHLhgCRQEp9GxbGjqc9/2RHJrExHpJI2WCVU8J/Bz1fWBYAeo2ht4HLpdRCTFqeXe2vrSj4O9iUbNiEiaUbi34hFGx0RaLyKSihTuIdydqtwTwm6rytE8MyKSPhTuIcyMglnfopa288ysP+Ur3Lvu9SRVJiLSOQr31sZ+idXDbqGysR/ugTszNWTl8czrB3TjDhFJGxotE0Efq6dpvrDjjlXxg94PkDPsTMzOSG5hIiJRUMs9hLtTWlbOp99cSl+rb7Et1+s0YkZE0obCPYSZkZ+Xw5CsyvA7VFWoW0ZE0oLCvZUbzh9BVU7kETOlq8pZvHZXgqsSEekchXsrZsbGomvD3pnpgbzLWbZxj06sikjK0wnVMGbM+yq+7WR8fSlU7eWYZ5HjdVxStYyi0/oye9YM3Z1JRFKaWu4RWPEcLDhTZC9rxAwKsw4ye9/d8Mpvkl2eiEi7FO7t8Ah3Z6pa9Q31u4tISlO4R+DuEe/ClF+/n+oa9buLSOpSuEdgZhzufWLYbQezB1NywWj1u4tIylK4tyN/5iK81d2Zjngugz//bQW7iKQ0hXt7gndn+ig7v3memVpyWf6SLmYSkdSmcG+Hu7P8pQrsaC1mYAbH24dMf/Mulj+8WAEvIilL4d4OM2Pquz9vM89MX6tn6rs/V9eMiKSsqMLdzKaZ2U4z221mt0TYZ46ZlZvZq2b269iWmTz5de9FXK+Wu4ikqg6vUDWzbGApMBWoAF40s5XuXh6yzwjgVmCiux8ys/CTs6SjgsLATbJb8YJCSsvKyc/LYcHUkUkoTEQksmha7ucAu939DXevBx4DLmq1z78BS939EIC7749tmUk0paTNiBnP6cPyAVexbIPmmRGR1BRNuA8BQpuuFcF1oUYCI81sg5k9b2bTYlVg0hXPwS5YghecigNHPQvqa/j0G0u557SdlMzSeHcRST3RTBwWLrlaN1V7ASOAyUAh8JyZjXH3D1o8kdk1wDUAQ4cO7XSxSVM8BwN81fX0Ck5HUGgHGbLvbuyVwsCQSRGRFBJNy70CODXkcSHwTph9fu/uDe7+JrCTQNi34O73u/sEd58wePDgrtacFJHmmXHdnUlEUlA04f4iMMLMhptZLnApsLLVPr8D/jeAmQ0i0E3zRiwLTab25pnR3ZlEJBV1GO7ufhSYD/wR2AE84e6vmlmpmV0Y3O2PQKWZlQN/Bha6e4R71aWf9uaZOdz7RPW5i0jKsWS1OidMmOCbN29Oymt3ybYn8FXXt+ia8Zw+2AVL1OcuIgljZlvcfUJH++kK1WgFR83Q5/jmVdarTzt/QUQkeRTuneRHQ06q1ryPr7oetj2RvIJERMJQuHdC9erbw46YqV59e5IqEhEJT+EeJXenf4R5ZvprnhkRSTEK9yiZWWCemXAKCjViRkRSisK9EyzCPDM2pUQtdxFJKQr3TvCxX2L5kK9T0TiIRg/OM9NQg6+/k+UPL2bx2l3JLlFEBFC4R83dKS0r52uvjeLX/a6gllx6WSMGWFUF09+8i6H7ytSCF5GUoHCPkpmRn5fDlROLWJjzeNi7M80+9KD63kUkJUQzK6QELZg6EnfH7twXdrtFmn9GRCTB1HLvJDPDI4ya8YLW09yLiCSHwr2T3J3lA67iiOe2WH/Ec/l+wyUsXrszSZWJiHxM4d5JZsZbp8zkJ8ddR2VjP5zAnUsasnpT8UEN1bVHdVJVRJJOfe5dcONnR9F44iiO/q4BC+Z4gR/mB70fIGfYmZidkdwCRaTHU8u9i7L+exG5XtdiXa7XYbozk4ikAIV7F3mEkTGR1ouIJJLCvQvcnarcE8Juq8o9QX3uIpJ0CvcuMDM2DruWeuvdYn299WbjsGt1IZOIJJ3CvYtmzPsqORf/BApODaywbHK8jhn7/1M37xCRpNNomW6w4L1Tm+6tagBVewOPQfdWFZGkUcu9m3R3JhFJRQr3btDdmUQkVSncu6H9uzMN0YlVEUkahXs32ZQSGrLyWqzznD4sH3C1bt4hIkmjcO+mxjFfZHHetYF5ZhzcofpoL57ZdYDq2gZ1zYhIUmi0TDe4O4tW76DygxqOy22gqROmwA9zd+4D5A4dr3lmRCQp1HLvhqa7M5X2+y15tJxnJo86sv57UZIqE5GeTuHeTQumjqSgfn/YbZpnRkSSReHeTZpnRkRSkcK9mzTPjIikoqjC3cymmdlOM9ttZre0s98XzczNbELsSkx9mmdGRFJNh+FuZtnAUmA6MBqYa2ajw+zXH7ge+J9YF5kOrHgOTCmBnD7gx5rnmWHV9Qp4EUm4aFru5wC73f0Nd68HHgMuCrPfIuBuoDaG9aWX9aXQap4ZGmoC60VEEiiacB8C7A15XBFc18zMzgROdfey9p7IzK4xs81mtvnAgQOdLjbV6e5MIpIqogn3cGcEm4eAmFkWsBj4WkdP5O73u/sEd58wePDg6KtMA+2NmtmfNYjFa3cmuCIR6cmiCfcK4NSQx4XAOyGP+wNjgKfNbA9wLrCyp51UNTPyZy6ilpajZhqBP9SPo7r2qIZFikjCRBPuLwIjzGy4meUClwIrmza6e5W7D3L3IncvAp4HLnT3zXGpOIVljbuE3AmX0RiS4VnA3JznKBn2qoZFikjCdBju7n4UmA/8EdgBPOHur5pZqZldGO8C04m7U72tjKxWGZ7rddj6UrXcRSRhopo4zN3XAGtarSuJsO/k7peVftyd0rJybq97L+xZCq+qoLSsnPy8HBZMHZn4AkWkR9EVqjFiZvTv3YsD2eFPFO/PGsSyDXs0DbCIJITCPYZu/OwotnzqujZTERzxXL5d+yWunFhEyazR6nsXkbhTuMdY6FQEDhz1LPpQz9d7PUHJsFeTXZ6I9BC6WUccWPEcHKhbPp88C8zzXmgHqV8xn1VbKni7cJb63UUkrtRyjwN3p2rVN9rcwCPX6/j0m0uprlG/u4jEl8I9DsyMgobwN/A4JauSkgvazLsmIhJTCvc4sYLCdteXlpWzeO2uRJYkIj2Iwj1eppTgOX1arGp0eKHXP1C6qlzDIkUkrhTuceJjv8QLBdNbTkVgMPZgGZXP/4orz9OwSBGJH42WiRMz4/TDG9tMRdDXAsMih1zwHQW7iMSNWu5xlF/3Xtj1p1glpWXl6pIRkbhRuMdThJOq1b1PYNmGPQp4EYkbhXs8Nd1TtZWC7AbuOW0n+Xk56poRkbhQn3s8Fc8J/HzqZrzm/ebJIq3mfWbvuxs7qxDQlaoiEntqucdb8RzIPa7NLMCmG2eLSBwp3BMh0g2ydeNsEYkThXsiRDixSp8Bia1DRHoMhXsiTCmBrJy26+s/hG1PJL4eEcl4CvdEKJ4Dvfu3XX+sXv3uIhIXCvdEqTkUfn1Vhca6i0jMKdwTJUK/uxcUaoZIEYk5hXuihJkl0oEjH1ZTuelXmiFSRGJKFzElSvCCJn/qZghe0GTAcceq+EHvB+g1dDxmZyS1RBHJHGq5J8jitbsofeuMsBc05Xod1atvV9eMiMSMwj0B3J3q2gaWbdgT8cKl/Lr96poRkZhRt0wCmBkls0bj7uzbPJDCrINt9jls/bh95umaSExEYkIt9wQxM755wRncfXQOdZ7dZvtxVkvW9ieTUJmIZCKFe4K4O6Vl5axs/Cc+ou00wL28AV9fqm4ZEYkJhXsCNAX7sg17+NfzhjHAPgq/Y1UFpas05l1Eui+qcDezaWa208x2m9ktYbbfaGblZrbNzNab2bDYl5q+zIz8vByunFjENy84g8O9Twy7XyNG5fMa8y4i3dfhCVUzywaWAlOBCuBFM1vp7uUhu/0VmODuR8zs/wJ3A5fEo+B0tWDqSNw9EPQzF+Grrg/M6R4im0Z+0PsBcoadqTHvItIt0bTczwF2u/sb7l4PPAZcFLqDu//Z3Y8EHz4PRJjjtmdrHglTPAe7YAlY2xOruV6HaTIxEemmaMJ9CLA35HFFcF0kVwNPdaeoHqF4Du6NYTd51d6w60VEohVNuIcbeB22Q9jMLgMmAN+PsP0aM9tsZpsPHDgQfZUZyN2pyj0h/DZgzaM/SmxBIpJRogn3CuDUkMeFwDutdzKz84HbgAvdvS7cE7n7/e4+wd0nDB48uCv1ZgwzY8OwawnXds8CzntrqU6qikiXRRPuLwIjzGy4meUClwIrQ3cwszOBnxMI9v2xLzMzzZz31bBfiwAK6j8+jAp5EemsDkfLuPtRM5sP/BHIBh5091fNrBTY7O4rCXTD9AN+Ezxp+La7XxjHujOGFZwK4frYLYvlDy/mrVNmcrjuKPl5OSyYOjLxBYpIWopqbhl3XwOsabWuJGT5/BjX1WP4lBIaVswnt1VPlvkxpr95Fz/Z/yHL3j+bKycWNQ+lFBHpiK5QTSJ3p/StM7ip7mqOhXkr+lo9//Lhw1w5sYiSWaMV7CISNYV7EjVduTrwHy8jK/wAJIbYQUqGvapgF5FO0ZS/SdZ05Sq7C8P2vZtB/Yr59HIna5wu+hWR6KjlniKWD7iKI54bdluu13Hg97exeO1OjZwRkago3FOAmfH2kFk8NfzWCJ0zMPjYAapqGjRrpIhEReGeIhZMHcnsKxYEhkaGYQZjK//Eso17NGukiHRI4Z5CzAymlBBuxgcDPv3mUo2cEZGoKNxTTfEcXCNnRKSbFO4pJjChWPibeZhBw4r5+LYnElyViKQbhXuKMTM2DruWeusddnuu13Hkd19Tn7uItEvhnoJmzPsqORf/JOLImb6N1Tz16yWAJhUTkfAU7inKiudApJEzwOdev4PGrY9RWqahkSLSlsI9hdmUkoit92waqf/d9VRu0g21RaQthXsqK56D9Tk+4uY86liU+wi3zzxdI2hEpAWFe4rz6d+LeHIVIJ/DZG1/MrCvWu8iEqSJw1JY05TAlXVXszj3PrLD3JTPAF/xFdydRW+P0U09RARQyz2lhU4JbJ//WcT+d/Nj1K+4Tv3vItJMLfcUt2DqSBobG1m0Gq7zfhxvH4bdr6n/vf/Mb6v/XUTUck8HWVlZ5Ofl8OfhX8Nz+kTcL5/D2PYnW7Tc1YoX6ZnUck8TgZt6LMBeKYQVXwE/1mYfA44t/3d+t6WCtwtnccP5IygtK1c/vEgPpJZ7GjEzKJ5DYzv979k0cvGeOznvte9QuqqcZRs0RbBIT6RwTzNNo2IOeb+I+2QZnHNwBddt+Sz3nLZTUwSL9EAK9zTTNIKmo/53MzjePmT2nlJY/bUEVigiqcCS9XV9woQJvnnz5qS8diZwd+yV3+ArvoKF6X9vsS9Q06uAvhf+AIrnJKZAEYkLM9vi7hM62k8t9zRlZvjYL7F82Ddo7OD3swF9j1bhy6/By24EPh5F0/qXu/rmRTKDwj2NNd1Y+8VBsyOeYG2xPw6bH+DIt4ay/OHFzLlvI6Vl5TQ2Bq58dXfNMimSITQUMs3dcP4ISmtu4tF3T+KuPr+i77HqMHdg/VhTK372m3cy22DfO4P4wY4ryB43h8N1R3lo41tcObGoRQteJ2NF0o/CPc2ZGfl9cuAfL6PvrO/A6q/hmx9oN+ADfy/wszDrIAs/ugc23sM+H0Sf46/gphnTmlvxW/dW8ZkRg1gwdWRzyLu7Al8kxemEaoZoeh9Ly8qp3PQr7sr7FX0b22/Fh38eaPpLh7wfdzRczusnTOecogEU9M2NeGGUmbUI/da/APQLQSQ2oj2hGlW4m9k04EdANvALd/9uq+29gV8CZwOVwCXuvqe951S4x8fitbuorm3g9pmn8+LSq/mHg8vJ6kamevA/TrC1H7oMOIbhVOeexLqT/523hszEDKprjnL7rMA8842NjXxrzWv0790fQDOrAAAHp0lEQVSLG84fEfYXQCovd/RYy3pPoqkX2jaCuiJm4W5m2cAuYCpQAbwIzHX38pB9/h9Q7O5fMbNLgYvd/ZL2nlfhHj+NjY2UlpXz0Ma3uDDrLyzKfYR8P0w3Pk9RCW31uwcXLfjLwEOWm/ZPm2XD3HFLhVq0HFhOt/ck0Ag63Psk1p7877w9ZFaXpwSJ5VDIc4Dd7v6Gu9cDjwEXtdrnIuDh4PKTwBTTd/CkycrKoqBPLuNPLWDguZfxwHnrKMud0eGQye6y4D80I3CVbNPjLFotk27Ljlmq1KLl9HxPHAPy6/7O9DfvYui+srgPO47mhOoQYG/I4wrg05H2cfejZlYFDAQOxqJI6bwFU0dyw/kjADAzFtsdrKj4B85/9+f0r/875sS9JS8ibfW1emYfehCzG+P6OtGEe7gIaP0rJ5p9MLNrgGsAhg4dGsVLS3eEfnlaMHUU7iMxu5HFa3cydN9qZh96AKoqQEEvklBWVRH314gm3CuAU0MeFwLvRNinwsx6AQXA+62fyN3vB+6HQJ97VwqWrmsK+5ZBv4uh+8o+DnqMpt/LynuR+PCCwrj/+4om3F8ERpjZcGAfcCnwL632WQlcAWwCvgj8t+s69pT2cdAH54kP+Yro7jz16yVMeuMH9D1WHVjHx2EfuqxWv0jnHPFcnhpwFbPjPDy4w3AP9qHPB/5IYCjkg+7+qpmVApvdfSXwAPCIme0m0GK/NG4VS8yFG641Y95Xcb/+422hV6wGf967bhdDK1Yz9e8/p1/d31tsS40RCrFaTreRGT1hOd3ek7ajZeI95kQXMUm3NI3Zbe9z5O5kZWXh7ikx5jjTx1T3hOV0e0+aJHKcu6YfkG5p+pC292EN3Sd0v3RZ7uixlhO/3NHjVFtur+540ayQIiIZSOEuIpKBFO4iIhlI4S4ikoEU7iIiGUjhLiKSgRTuIiIZKGkXMZnZAeCtbj7NIFJv5slUrAlUV2ekYk2gujorFeuKRU3D3H1wRzslLdxjwcw2R3OlViKlYk2gujojFWsC1dVZqVhXImtSt4yISAZSuIuIZKB0D/f7k11AGKlYE6iuzkjFmkB1dVYq1pWwmtK6z11ERMJL95a7iIiEkZbhbmbTzGynme02s1uSWMepZvZnM9thZq+a2VeD6+8ws31mtjX4Z0YSattjZq8EX39zcN3xZrbWzF4P/hyQwHpGhRyPrWZWbWY3JONYmdmDZrbfzLaHrAt7bCxgSfCzts3MzkpwXd83s9eCr73CzD4RXF9kZjUhx+2+BNYU8T0zs1uDx2qnmX0uHjW1U9fjITXtMbOtwfUJOVbB14qUCYn/fDXdQCFd/hC4G9TfgE8CucDLwOgk1XIycFZwuT+wCxgN3AHclOTjtAcY1Grd3cAtweVbgO8l8T38OzAsGccK+AxwFrC9o2MDzACeInBTnXOB/0lwXZ8FegWXvxdSV1HofgmuKex7Fvzsvwz0BoYH/51mJ6quVtvvAUoSeayCrxUpExL++UrHlvs5wG53f8Pd64HHgIuSUYi7v+vuLwWXDwM7gCHJqCVKFwEPB5cfBj6fpDqmAH9z9+5exNYl7v4sbW/gHunYXAT80gOeBz5hZicnqi53/5O7Hw0+fJ7ADeoTJsKxiuQi4DF3r3P3N4HdBP69JrQuC9wNYw7wX/F47fa0kwkJ/3ylY7gPAfaGPK4gBQLVzIqAM4H/Ca6aH/ya9WAiuz9COPAnM9tiZtcE153o7u9C4EMInJCEuiBwj93Qf3jJPlYQ+dik0uftKgKtvCbDzeyvZvaMmf1zgmsJ956lyrH6Z+A9d389ZF3Cj1WrTEj45ysdwz3cPaqSOuTHzPoBvwVucPdq4GfAp4DxwLsEviIm2kR3PwuYDlxrZp9JQg1tmFkucCHwm+CqVDhW7UmJz5uZ3QYcBR4NrnoXGOruZwI3Ar82s/wElRPpPUuJYwXMpWXjIeHHKkwmRNw1zLqYHLN0DPcK4NSQx4XAO0mqBTPLIfAmPuruywHc/T13P+bujcB/Eqevpu1x93eCP/cDK4I1vNf0lS/4c3+i6yLwy+Yld38vWF/Sj1VQpGOT9M+bmV0BzALmebCjNtj1URlc3kKgf3tkIupp5z1LhWPVC5gNPN60LtHHKlwmkITPVzqG+4vACDMbHmwFXgqsTEYhwb69B4Ad7v7DkPWhfWYXA9tb/90413WcmfVvWiZwUm47geN0RXC3K4DfJ7KuoBatqmQfqxCRjs1K4PLgqIZzgaqmr9eJYGbTgJuBC939SMj6wWaWHVz+JDACeCNBNUV6z1YCl5pZbzMbHqzphUTUFOJ84DV3r2hakchjFSkTSMbnKxFnkGP9h8AZ5l0EfgPflsQ6/onAV6htwNbgnxnAI8ArwfUrgZMTXNcnCYxaeBl4tekYAQOB9cDrwZ/HJ7iuvkAlUBCyLuHHisAvl3eBBgItp6sjHRsCX5uXBj9rrwATElzXbgJ9sk2fr/uC+34h+N6+DLwEXJDAmiK+Z8BtwWO1E5ieyGMVXP8Q8JVW+ybkWAVfK1ImJPzzpStURUQyUDp2y4iISAcU7iIiGUjhLiKSgRTuIiIZSOEuIpKBFO4iIhlI4S4ikoEU7iIiGej/A8Psju6vAa8xAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2a3539d8828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_people_set=[]\n",
    "prob_set=[]\n",
    "prob_set_theory=[]\n",
    "for i in range(1,200):\n",
    "    num_people_set.append(i)\n",
    "    prob_set.append(find_prob(1000,i))\n",
    "    prob_set_theory.append(P2(i))\n",
    "\n",
    "plt.scatter(num_people_set,prob_set,marker='x')\n",
    "plt.scatter(num_people_set,prob_set_theory,marker='o')\n",
    "plt.legend([\"simulation\",\"theory\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# problem3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.misc import factorial\n",
    "def nCr(n,r):\n",
    "    f = factorial\n",
    "    return f(n) / f(r) / f(n-r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Anaconda\\lib\\site-packages\\ipykernel_launcher.py:4: DeprecationWarning: `factorial` is deprecated!\n",
      "Importing `factorial` from scipy.misc is deprecated in scipy 1.0.0. Use `scipy.special.factorial` instead.\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.7272727272727273"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nCr(8,4)*nCr(4,2)/nCr(12,6)+nCr(8,5)*4/nCr(12,6)+nCr(8,6)/nCr(12,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(simulation)The probability that cards include exactly 3 Kings is 0.006100\n",
      "(mathematically)The probability that cards include exactly 3 Kings is 0.005818\n"
     ]
    }
   ],
   "source": [
    "card=np.arange(52)\n",
    "king_set=np.array([0,1,2,3])\n",
    "count=0\n",
    "n=10000\n",
    "for i in range(n):\n",
    "    np.random.shuffle(card)\n",
    "    if len(np.intersect1d(king_set,card[0:7]))==3:\n",
    "        count=count+1\n",
    "print(\"(simulation)The probability that cards include exactly 3 Kings is %f\" %(count/n))\n",
    "print(\"(mathematically)The probability that cards include exactly 3 Kings is %f\" %(nCr(4,3)*nCr(48,4)/nCr(52,7)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mathematical result: $ area = ln2 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.684"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point_set=np.random.uniform(0,1,[1000,2])\n",
    "np.sum(1/(1+point_set[:,0])>point_set[:,1])/1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mathematical soluntion is \n",
      "p1= 0.5 \n",
      "p2= 0.25 \n",
      "p3= 0.5 \n",
      "p4= 0 \n",
      "p5= 0.5\n",
      "The simulation soluntion is \n",
      "p1= 0.509808 \n",
      "p2= 0.242479 \n",
      "p3= 0.509835 \n",
      "p4= 0.000000 \n",
      "p5= 0.506143\n"
     ]
    }
   ],
   "source": [
    "n=10000\n",
    "toss_set=pd.DataFrame(np.floor(np.random.uniform(0,2,[n,3])),columns=list('ABC'))\n",
    "p1=(toss_set[(toss_set.A==1)].sum(axis=1)==2).sum()/np.size(toss_set[(toss_set.A==1)],0)\n",
    "p2=(toss_set[(toss_set.A==0)].sum(axis=1)==2).sum()/np.size(toss_set[(toss_set.A==0)],0)\n",
    "p3=(toss_set[(toss_set.A==1)&(toss_set.B==1)].sum(axis=1)==2).sum()/np.size(toss_set[(toss_set.A==1)&(toss_set.B==1)],0)\n",
    "p4=(toss_set[(toss_set.A==0)&(toss_set.B==0)].sum(axis=1)==2).sum()/np.size(toss_set[(toss_set.A==0)&(toss_set.B==0)],0)\n",
    "p5=(toss_set[(toss_set.A==1)&(toss_set.C==1)].sum(axis=1)==2).sum()/np.size(toss_set[(toss_set.A==1)&(toss_set.C==1)],0)\n",
    "print(\"The mathematical soluntion is \\np1= 0.5 \\np2= 0.25 \\np3= 0.5 \\np4= 0 \\np5= 0.5\")\n",
    "print(\"The simulation soluntion is \\np1= %f \\np2= %f \\np3= %f \\np4= %f \\np5= %f\" %(p1,p2,p3,p4,p5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mathematical soluntion is p= 0.533333\n",
      "The simulation soluntion is p= 0.539100\n"
     ]
    }
   ],
   "source": [
    "print(\"The mathematical soluntion is p= %f\" %(nCr(4,1)*nCr(2,1)/nCr(6,2)))\n",
    "count=0\n",
    "n=10000\n",
    "a=np.hstack((np.zeros(4),np.ones(2)))\n",
    "for i in range(n):\n",
    "    np.random.shuffle(a)\n",
    "    if a[:2].sum()==1:\n",
    "        count=count+1\n",
    "print(\"The simulation soluntion is p= %f\" %(count/n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "when N= 1, log(N!)= 0.000000 its approximation is 0.324404\n",
      "when N= 2, log(N!)= 0.693147 its approximation is 0.751520\n",
      "when N= 3, log(N!)= 1.791759 its approximation is 1.677227\n",
      "when N= 4, log(N!)= 3.178054 its approximation is 2.935248\n",
      "when N= 5, log(N!)= 4.787492 its approximation is 4.442679\n",
      "when N= 6, log(N!)= 6.579251 its approximation is 6.149752\n",
      "when N= 7, log(N!)= 8.525161 its approximation is 8.023260\n",
      "when N= 8, log(N!)= 10.604603 its approximation is 10.039468\n",
      "when N= 9, log(N!)= 12.801827 its approximation is 12.180565\n",
      "when N= 10, log(N!)= 15.104413 its approximation is 14.432691\n",
      "when N= 11, log(N!)= 17.502308 its approximation is 16.784756\n",
      "when N= 12, log(N!)= 19.987214 its approximation is 19.227682\n",
      "when N= 13, log(N!)= 22.552164 its approximation is 21.753904\n",
      "when N= 14, log(N!)= 25.191221 its approximation is 24.357020\n",
      "when N= 15, log(N!)= 27.899271 its approximation is 27.031539\n",
      "when N= 16, log(N!)= 30.671860 its approximation is 29.772705\n",
      "when N= 17, log(N!)= 33.505073 its approximation is 32.576354\n",
      "when N= 18, log(N!)= 36.395445 its approximation is 35.438812\n",
      "when N= 19, log(N!)= 39.339884 its approximation is 38.356813\n",
      "when N= 20, log(N!)= 42.335616 its approximation is 41.327436\n",
      "when N= 21, log(N!)= 45.380139 its approximation is 44.348050\n",
      "when N= 22, log(N!)= 48.471181 its approximation is 47.416275\n",
      "when N= 23, log(N!)= 51.606676 its approximation is 50.529948\n",
      "when N= 24, log(N!)= 54.784729 its approximation is 53.687093\n",
      "when N= 25, log(N!)= 58.003605 its approximation is 56.885900\n",
      "when N= 26, log(N!)= 61.261702 its approximation is 60.124702\n",
      "when N= 27, log(N!)= 64.557539 its approximation is 63.401961\n",
      "when N= 28, log(N!)= 67.889743 its approximation is 66.716253\n",
      "when N= 29, log(N!)= 71.257039 its approximation is 70.066256\n",
      "when N= 30, log(N!)= 74.658236 its approximation is 73.450739\n",
      "when N= 31, log(N!)= 78.092224 its approximation is 76.868552\n",
      "when N= 32, log(N!)= 81.557959 its approximation is 80.318621\n",
      "when N= 33, log(N!)= 85.054467 its approximation is 83.799938\n",
      "when N= 34, log(N!)= 88.580828 its approximation is 87.311556\n",
      "when N= 35, log(N!)= 92.136176 its approximation is 90.852583\n",
      "when N= 36, log(N!)= 95.719695 its approximation is 94.422180\n",
      "when N= 37, log(N!)= 99.330612 its approximation is 98.019553\n",
      "when N= 38, log(N!)= 102.968199 its approximation is 101.643952\n",
      "when N= 39, log(N!)= 106.631760 its approximation is 105.294665\n",
      "when N= 40, log(N!)= 110.320640 its approximation is 108.971017\n",
      "when N= 41, log(N!)= 114.034212 its approximation is 112.672369\n",
      "when N= 42, log(N!)= 117.771881 its approximation is 116.398110\n",
      "when N= 43, log(N!)= 121.533082 its approximation is 120.147659\n",
      "when N= 44, log(N!)= 125.317271 its approximation is 123.920463\n",
      "when N= 45, log(N!)= 129.123934 its approximation is 127.715993\n",
      "when N= 46, log(N!)= 132.952575 its approximation is 131.533745\n",
      "when N= 47, log(N!)= 136.802723 its approximation is 135.373235\n",
      "when N= 48, log(N!)= 140.673924 its approximation is 139.234001\n",
      "when N= 49, log(N!)= 144.565744 its approximation is 143.115599\n",
      "when N= 50, log(N!)= 148.477767 its approximation is 147.017605\n",
      "when N= 51, log(N!)= 152.409593 its approximation is 150.939611\n",
      "when N= 52, log(N!)= 156.360836 its approximation is 154.881223\n",
      "when N= 53, log(N!)= 160.331128 its approximation is 158.842066\n",
      "when N= 54, log(N!)= 164.320112 its approximation is 162.821776\n",
      "when N= 55, log(N!)= 168.327445 its approximation is 166.820005\n",
      "when N= 56, log(N!)= 172.352797 its approximation is 170.836414\n",
      "when N= 57, log(N!)= 176.395848 its approximation is 174.870681\n",
      "when N= 58, log(N!)= 180.456291 its approximation is 178.922490\n",
      "when N= 59, log(N!)= 184.533829 its approximation is 182.991541\n",
      "when N= 60, log(N!)= 188.628173 its approximation is 187.077540\n",
      "when N= 61, log(N!)= 192.739047 its approximation is 191.180206\n",
      "when N= 62, log(N!)= 196.866182 its approximation is 195.299265\n",
      "when N= 63, log(N!)= 201.009316 its approximation is 199.434453\n",
      "when N= 64, log(N!)= 205.168199 its approximation is 203.585513\n",
      "when N= 65, log(N!)= 209.342587 its approximation is 207.752198\n",
      "when N= 66, log(N!)= 213.532241 its approximation is 211.934267\n",
      "when N= 67, log(N!)= 217.736934 its approximation is 216.131488\n",
      "when N= 68, log(N!)= 221.956442 its approximation is 220.343633\n",
      "when N= 69, log(N!)= 226.190548 its approximation is 224.570484\n",
      "when N= 70, log(N!)= 230.439044 its approximation is 228.811828\n",
      "when N= 71, log(N!)= 234.701723 its approximation is 233.067457\n",
      "when N= 72, log(N!)= 238.978390 its approximation is 237.337171\n",
      "when N= 73, log(N!)= 243.268849 its approximation is 241.620773\n",
      "when N= 74, log(N!)= 247.572914 its approximation is 245.918074\n",
      "when N= 75, log(N!)= 251.890402 its approximation is 250.228888\n",
      "when N= 76, log(N!)= 256.221136 its approximation is 254.553035\n",
      "when N= 77, log(N!)= 260.564941 its approximation is 258.890340\n",
      "when N= 78, log(N!)= 264.921650 its approximation is 263.240631\n",
      "when N= 79, log(N!)= 269.291098 its approximation is 267.603743\n",
      "when N= 80, log(N!)= 273.673124 its approximation is 271.979513\n",
      "when N= 81, log(N!)= 278.067573 its approximation is 276.367783\n",
      "when N= 82, log(N!)= 282.474293 its approximation is 280.768399\n",
      "when N= 83, log(N!)= 286.893133 its approximation is 285.181209\n",
      "when N= 84, log(N!)= 291.323950 its approximation is 289.606067\n",
      "when N= 85, log(N!)= 295.766601 its approximation is 294.042830\n",
      "when N= 86, log(N!)= 300.220949 its approximation is 298.491358\n",
      "when N= 87, log(N!)= 304.686857 its approximation is 302.951514\n",
      "when N= 88, log(N!)= 309.164194 its approximation is 307.423163\n",
      "when N= 89, log(N!)= 313.652830 its approximation is 311.906176\n",
      "when N= 90, log(N!)= 318.152640 its approximation is 316.400425\n",
      "when N= 91, log(N!)= 322.663499 its approximation is 320.905785\n",
      "when N= 92, log(N!)= 327.185288 its approximation is 325.422134\n",
      "when N= 93, log(N!)= 331.717887 its approximation is 329.949352\n",
      "when N= 94, log(N!)= 336.261182 its approximation is 334.487323\n",
      "when N= 95, log(N!)= 340.815059 its approximation is 339.035932\n",
      "when N= 96, log(N!)= 345.379407 its approximation is 343.595067\n",
      "when N= 97, log(N!)= 349.954118 its approximation is 348.164619\n",
      "when N= 98, log(N!)= 354.539086 its approximation is 352.744480\n",
      "when N= 99, log(N!)= 359.134205 its approximation is 357.334545\n",
      "when N= 100, log(N!)= 363.739376 its approximation is 361.934711\n",
      "The exact value of log(N!) for N=200 is 863.231987\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,101):\n",
    "    i=float(i)\n",
    "    print(\"when N= %d, log(N!)= %f its approximation is %f\" %(i,np.log(float(math.factorial(i))),i*np.log(i+0.5)+np.log(np.sqrt(2*np.pi)*pow(math.e,-i))))\n",
    "print(\"The exact value of log(N!) for N=200 is %f\" %(np.log(np.arange(1,201)).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nCr(n,r):\n",
    "    return np.exp(np.log(np.arange(1,n+1)).sum()-np.log(np.arange(1,r+1)).sum()-np.log(np.arange(1,n-r+1)).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "When k= 95, P_bino(k)= 0.163918  P_hype(k)= 0.184761\n",
      "When k= 96, P_bino(k)= 0.133752  P_hype(k)= 0.142656\n",
      "When k= 97, P_bino(k)= 0.086410  P_hype(k)= 0.085604\n",
      "When k= 98, P_bino(k)= 0.041442  P_hype(k)= 0.037442\n",
      "When k= 99, P_bino(k)= 0.013116  P_hype(k)= 0.010615\n",
      "When k= 100, P_bino(k)= 0.002055  P_hype(k)= 0.001463\n"
     ]
    }
   ],
   "source": [
    "N=1000\n",
    "M=100\n",
    "p=0.94\n",
    "k=[]\n",
    "bino=[]\n",
    "hype=[]\n",
    "for i in range(95,101):\n",
    "    k.append(i)\n",
    "    bino.append(nCr(M,i)*np.power(p,i)*np.power(1-p,M-i))\n",
    "    hype.append(nCr(p*N,i)*nCr((1-p)*N,M-i)/nCr(N,M))\n",
    "    print(\"When k= %d, P_bino(k)= %f  P_hype(k)= %f\" %(i,nCr(M,i)*np.power(p,i)*np.power(1-p,M-i),nCr(p*N,i)*nCr((1-p)*N,M-i)/nCr(N,M)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# problem9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mathematical resulf of p(x>=4) = 0.421875\n",
      "The simulation resulf of p(x>=4) = 0.416000\n"
     ]
    }
   ],
   "source": [
    "print(\"The mathematical resulf of p(x>=4) = %f\" %(np.power(0.75,3)))\n",
    "n=1000\n",
    "sample=np.random.geometric(0.25,n)\n",
    "print(\"The simulation resulf of p(x>=4) = %f\" % (sum(sample>3)/n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# problem10(f,c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Anaconda\\lib\\site-packages\\ipykernel_launcher.py:5: DeprecationWarning: `factorial` is deprecated!\n",
      "Importing `factorial` from scipy.misc is deprecated in scipy 1.0.0. Use `scipy.special.factorial` instead.\n",
      "  \"\"\"\n",
      "E:\\Anaconda\\lib\\site-packages\\ipykernel_launcher.py:4: DeprecationWarning: `factorial` is deprecated!\n",
      "Importing `factorial` from scipy.misc is deprecated in scipy 1.0.0. Use `scipy.special.factorial` instead.\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c88476aeb8>"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGyBJREFUeJzt3X9wX3Wd7/Hnq2lDiwiUNlxvm8SGvQWpECl8raxwkbFoy+q07AhYFG1VbkSpuHoXbxkc9BZndO0OqDtFqAiXRaVUlG5WVztukfUyA9j0B/3FZgm1tmnU/pK61xbalPf943vKfPMlISfJN/m257weM5l8z+d8Pt/v+8zJvHLyOSfnKCIwM7N8GFXtAszMbOQ49M3McsShb2aWIw59M7McceibmeWIQ9/MLEdShb6k2ZLaJXVIWvQ6/a6WFJIKJW23JuPaJc2qRNFmZjY4o/vrIKkGWAq8B+gE1khqjYitZf3eCNwMPFPSNg2YB7wVmAT8q6SzI+Jo5TbBzMzSSnOkPwPoiIhtEXEYWA7M7aXfHcDXgZdK2uYCyyPi5Yj4DdCRvJ+ZmVVBv0f6wGRgZ8lyJ/CO0g6SpgMNEfETSX9bNvbpsrGTX+/DJk6cGFOmTElRlpmZHbN27dq9EVHXX780oa9e2l69d4OkUcBdwIKBji15jxagBaCxsZG2trYUZZmZ2TGSfpumX5rpnU6goWS5HugqWX4jcB7whKTtwMVAa3Iyt7+xAETEsogoREShrq7fX1RmZjZIaUJ/DTBVUpOkWoonZluPrYyIAxExMSKmRMQUitM5cyKiLek3T9JJkpqAqcCvK74VZmaWSr/TOxHRLWkhsAqoAe6PiC2SFgNtEdH6OmO3SFoBbAW6gZt85Y6ZWfXoeLu1cqFQCM/pm9lgHTlyhM7OTl566aX+O5+Axo4dS319PWPGjOnRLmltRBT6GPaqNCdyzcxOGJ2dnbzxjW9kypQpSL1dS3Liigj27dtHZ2cnTU1Ng3oP34bBzDLlpZdeYsKECZkLfABJTJgwYUh/xTj0zSxzshj4xwx12xz6ZmY54tA3M6uyG264ga1bt/bfsQJ8ItfMrMruu+++EfssH+mbWa6tXL+LS772OE2LfsolX3uclet3Dfk9t2/fzlve8hbmz59Pc3MzV199NQcPHmT16tVMnz6d888/n49//OO8/PLLAFx++eW0tbVx9OhRFixYwHnnncf555/PXXfdNeRaymUz9DeugLvOgy+fXvy+cUW1KzKz49DK9bu49ceb2PXiIQLY9eIhbv3xpooEf3t7Oy0tLWzcuJFTTz2VO++8kwULFvDII4+wadMmuru7+fa3v91jzIYNG9i1axebN29m06ZNfOxjHxtyHeWyF/obV9D9T5+BAzuBgAM7i8sOfjMrs2RVO4eO9LxJwKEjR1myqn3I793Q0MAll1wCwPXXX8/q1atpamri7LPPBmD+/Pn86le/6jHmrLPOYtu2bXzmM5/h5z//OaeeeuqQ6yiXudA/+LPbGX205zWso4++xMGf3V6liszseNX14qEBtQ/EYC6tHD9+PM8++yyXX345S5cu5YYbbhhyHeUyF/pjD/1+QO1mll+TTh83oPaB2LFjB0899RQADz/8MFdccQXbt2+no6MDgIceeoh3vetdPcbs3buXV155hQ984APccccdrFu3bsh1lMtc6He9MmFA7WaWX7fMOodxY2p6tI0bU8Mts84Z8nufe+65PPjggzQ3N7N//34+97nP8cADD3DNNddw/vnnM2rUKG688cYeY3bt2sXll1/OBRdcwIIFC/jqV7865DrKZe6Szftqr+cLR+7mZB1+te1g1HJf7fV8uXplmdlx6KrpxQf5LVnVTteLh5h0+jhumXXOq+1DMWrUKO65554ebTNnzmT9+vWv6fvEE0+8+no4ju5LZS70L3hfC7c/1s3fxHImaR9dMYFvMI9L39dS7dLM7Dh01fTJFQn5E0XmQr+48z7NB1fNrPhvbjOzNKZMmcLmzZurXUavMhf6kL/f3GZmaaU6kStptqR2SR2SFvWy/kZJmyRtkPSkpGlJ+xRJh5L2DZLuee27m5nZSOn3SF9SDbAUeA/FB52vkdQaEaV3B/pBRNyT9J8D3AnMTta9EBEXVLZsMzMbjDRH+jOAjojYFhGHgeXA3NIOEfGnksU3AMfXMxjNzAxIF/qTgZ0ly51JWw+SbpL0AvB14OaSVU2S1kv6N0n/fUjVmpmdALZv38555533mvaRvIVyX9KcyO3tf4lfcyQfEUuBpZI+BHwRmA/8DmiMiH2SLgJWSnpr2V8GSGoBWgAaGxsHuAlmZieGkbyFcl/SHOl3Ag0ly/VA1+v0Xw5cBRARL0fEvuT1WuAF4OzyARGxLCIKEVGoq6tLW7uZ2dAN0115u7u7X3Nr5WO3UAY45ZRTuO2223jb297GxRdfzB/+8AcAfvvb3zJz5kyam5uZOXMmO3bsqEg9x6QJ/TXAVElNkmqBeUBraQdJU0sW3wc8n7TXJSeCkXQWMBXYVonCzcyGbOMK+Oebe9yVl3++uSLBX35r5bvvvrvH+j//+c9cfPHFPPvss1x22WV85zvfAWDhwoV89KMfZePGjXz4wx/m5ptv7u3tB63f0I+IbmAhsAp4DlgREVskLU6u1AFYKGmLpA3A5ylO7QBcBmyU9CzwKHBjROyv6BaYmQ3W6sVwpOyOmkcOFduHqPzWyk8++WSP9bW1tbz//e8H4KKLLmL79u0APPXUU3zoQx8C4CMf+chrxg1Vqn/Oioh/Af6lrO32ktef7WPcj4AfDaVAM7Nhc6BzYO0DUH5r5fLlMWPGvNpWU1NDd3d3qvcZqszdZdPMLLXT6gfWPgDlt1a+9NJLU4175zvfyfLlywH4/ve/n3pcWg59M8uvmbfDmLJ7548ZV2wfovJbK3/qU59KNe5b3/oWDzzwAM3NzTz00EN885vfHHItpRRxfP0fVaFQiGNnt83MBuq5557j3HPPTT9g44riHP6BzuIR/szbofna4SuwAnrbRklrI6LQ39hM3nDNzCy15muP+5CvJE/vmJnliEPfzDLneJu2rqShbptD38wyZezYsezbty+TwR8R7Nu3j7Fjxw76PTynb2aZUl9fT2dnJ3v27Kl2KcNi7Nix1NcP/pJSh76ZZcqYMWNoamqqdhnHLU/vmJnliEPfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjqUJf0mxJ7ZI6JC3qZf2NkjZJ2iDpSUnTStbdmoxrlzSrksWbmdnA9Bv6yTNulwJXAtOA60pDPfGDiDg/Ii4Avg7cmYydRvGZum8FZgN3H3tmrpmZjbw0R/ozgI6I2BYRh4HlwNzSDhHxp5LFNwDHbnoxF1geES9HxG+AjuT9zMysCtLchmEysLNkuRN4R3knSTdRfCh6LfDukrFPl42dPKhKzcxsyNIc6ff2VN7X3L4uIpZGxF8A/wv44kDGSmqR1CapLas3STIzOx6kCf1OoKFkuR7oep3+y4GrBjI2IpZFRCEiCnV1dSlKMjOzwUgT+muAqZKaJNVSPDHbWtpB0tSSxfcBzyevW4F5kk6S1ARMBX499LLNzGww+p3Tj4huSQuBVUANcH9EbJG0GGiLiFZgoaQrgCPAH4H5ydgtklYAW4Fu4KaIODpM22JmZv3Q8fZ0mUKhEG1tbdUuw8zshCJpbUQU+uvn/8g1M8sRh76ZWY449M3McsShb2aWIw59M7McceibmeWIQ9/MLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliMOfTOzHHHom5nliEPfzCxHHPpmZjmSKvQlzZbULqlD0qJe1n9e0lZJGyWtlvTmknVHJW1IvlrLx5qZ2cjp93GJkmqApcB7KD7ofI2k1ojYWtJtPVCIiIOSPgV8Hfhgsu5QRFxQ4brNzGwQ0hzpzwA6ImJbRBwGlgNzSztExC8j4mCy+DRQX9kyzcysEtKE/mRgZ8lyZ9LWl08APytZHiupTdLTkq4aRI1mZlYh/U7vAOqlrdenqUu6HigA7yppboyILklnAY9L2hQRL5SNawFaABobG1MVbmZmA5fmSL8TaChZrge6yjtJugK4DZgTES8fa4+IruT7NuAJYHr52IhYFhGFiCjU1dUNaAPMzCy9NKG/BpgqqUlSLTAP6HEVjqTpwL0UA393Sft4SSclrycClwClJ4DNzGwE9Tu9ExHdkhYCq4Aa4P6I2CJpMdAWEa3AEuAU4IeSAHZExBzgXOBeSa9Q/AXztbKrfszMbAQpotfp+aopFArR1tZW7TLMzE4oktZGRKG/fv6PXDOzHHHom5nliEPfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjDn0zsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY449M3McsShb2aWI6lCX9JsSe2SOiQt6mX95yVtlbRR0mpJby5ZN1/S88nX/EoWb2ZmA9Nv6EuqAZYCVwLTgOskTSvrth4oREQz8Cjw9WTsGcCXgHcAM4AvSRpfufLNzGwg0hzpzwA6ImJbRBwGlgNzSztExC8j4mCy+DRQn7yeBfwiIvZHxB+BXwCzK1O6mZkNVJrQnwzsLFnuTNr68gngZ4Mca2Zmw2h0ij7qpa3Xp6lLuh4oAO8ayFhJLUALQGNjY4qSzMxsMNIc6XcCDSXL9UBXeSdJVwC3AXMi4uWBjI2IZRFRiIhCXV1d2trNzGyA0oT+GmCqpCZJtcA8oLW0g6TpwL0UA393yapVwHsljU9O4L43aTMzsyrod3onIrolLaQY1jXA/RGxRdJioC0iWoElwCnADyUB7IiIORGxX9IdFH9xACyOiP3DsiVmZtYvRfQ6PV81hUIh2traql2GmdkJRdLaiCj01y/NidwT2prWe2lYt4QzYw+7VcfOC2/h7XM+We2yzMyqItOhv6b1Xs5b+0XG6TAI3sQeTlv7RdaAg9/McinT995pWLekGPglxukwDeuWVKkiM7PqynTonxl7+mjfO8KVmJkdHzId+rvV+zX/uzVxhCsxMzs+ZDr0d154C4eitkfboahl54W3VKkiM7PqynTov33OJ9l80Vf4PXW8EuL31LH5oq/4JK6Z5Zav0zczy4C01+ln+kjfzMx6cuibmeWIQ9/MLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliOpQl/SbEntkjokLepl/WWS1knqlnR12bqjkjYkX63lY83MbOT0ez99STXAUuA9FB90vkZSa0RsLem2A1gA/G0vb3EoIi6oQK1mZjZEaR6iMgPoiIhtAJKWA3OBV0M/IrYn614ZhhrNzKxC0kzvTAZ2lix3Jm1pjZXUJulpSVcNqDozM6uoNEf66qVtIHdpa4yILklnAY9L2hQRL/T4AKkFaAFobGwcwFubmdlApDnS7wQaSpbrga60HxARXcn3bcATwPRe+iyLiEJEFOrqen/wiZmZDV2a0F8DTJXUJKkWmAekugpH0nhJJyWvJwKXUHIuwMzMRla/oR8R3cBCYBXwHLAiIrZIWixpDoCkt0vqBK4B7pW0JRl+LtAm6Vngl8DXyq76MTOzEeSHqJiZZYAfomJmZq/h0DczyxGHvplZjjj0zcxyxKFvZpYjDn0zsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY449M3McsShb2aWIw59M7McceibmeWIQ9/MLEdShb6k2ZLaJXVIWtTL+sskrZPULenqsnXzJT2ffM2vVOFmZjZw/Ya+pBpgKXAlMA24TtK0sm47gAXAD8rGngF8CXgHMAP4kqTxQy/bzMwGI82R/gygIyK2RcRhYDkwt7RDRGyPiI3AK2VjZwG/iIj9EfFH4BfA7ArUbWZmg5Am9CcDO0uWO5O2NIYy1szMKixN6KuXtrRPU081VlKLpDZJbXv27En51mZmNlBpQr8TaChZrge6Ur5/qrERsSwiChFRqKurS/nWZmY2UGlCfw0wVVKTpFpgHtCa8v1XAe+VND45gfvepM3MzKqg39CPiG5gIcWwfg5YERFbJC2WNAdA0tsldQLXAPdK2pKM3Q/cQfEXxxpgcdJmZmZVoIi00/Mjo1AoRFtbW7XLMDM7oUhaGxGF/vr5P3LNzHLEoW9mliMOfTOzHHHom5nliEPfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjDn0zsxwZXe0CRtLK9btYsqqdrhcPMen0cdwy6xyumu5nuphZfuQm9Feu38WTj93NIyxn0kl76To4kW88Ng/4tIPfzHIjN9M7G366jMVaRv2ovYwS1I/ay2ItY8NPl1W7NDOzEZOb0L/h8Pc4WYd7tJ2sw9xw+HtVqsjMbOTlJvQnjdo3oHYzsyxKFfqSZktql9QhaVEv60+S9Eiy/hlJU5L2KZIOSdqQfN1T2fLTe2ncmwbUbmaWRf2GvqQaYClwJTANuE7StLJunwD+GBH/DbgL+LuSdS9ExAXJ140VqnvATr5yMd01Y3u0ddeM5eQrF1epIjOzkZfmSH8G0BER2yLiMLAcmFvWZy7wYPL6UWCmJFWuzApovpbRc/8BTmsABKc1FJebr612ZWZmIybNJZuTgZ0ly53AO/rqExHdkg4AE5J1TZLWA38CvhgR/3doJQ9B87UOeTPLtTSh39sRe/nT1Pvq8zugMSL2SboIWCnprRHxpx6DpRagBaCxsTFFSWZmNhhppnc6gYaS5Xqgq68+kkYDpwH7I+LliNgHEBFrgReAs8s/ICKWRUQhIgp1dXUD3wozM0slTeivAaZKapJUC8wDWsv6tALzk9dXA49HREiqS04EI+ksYCqwrTKlm5nZQPU7vZPM0S8EVgE1wP0RsUXSYqAtIlqB7wIPSeoA9lP8xQBwGbBYUjdwFLgxIvYPx4aYmVn/FFE+PV9dhUIh2traql2GmdkJRdLaiCj01y83/5FrZmYOfTOzXHHom5nliEPfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjae6ymUkr1+9iyap2ul48xKTTx3HLrHO4avrkapdlZjaschn6K9fv4snH7uYRljPppL10HZzINx6bB3zawW9mmZbL6Z0NP13GYi2jftReRgnqR+1lsZax4afLql2amdmwymXo33D4e5yswz3aTtZhbjj8vSpVZGY2MnIZ+pNG7RtQu5lZVuQy9F8a96YBtZuZZUUuQ//kKxfTXTO2R9sharn1wF9zydceZ+X6XVWqzMxseKUKfUmzJbVL6pC0qJf1J0l6JFn/jKQpJetuTdrbJc2qXOlD0Hwto+f+A5zWQCD2xykcilruGnM3jxz8Hzz52N0OfjPLpH6fnJU84/Y/gPdQfAD6GuC6iNha0ufTQHNE3ChpHvDXEfFBSdOAh4EZwCTgX4GzI+JoX5830k/O+vJXvsQXjtzd48TuKwECXtQpgDgt/pMDfu3XVXq9W3X85oxLadr/JGfGnuOiJr+u/D7eeeEtvH3OJ183r15P2idnpQn9vwS+HBGzkuVbASLiqyV9ViV9npI0Gvg9UAcsKu1b2q+vzxvp0O+8/S+oH7V3xD7PbDAiQKp2FTacDkUtmy/6yqCDv5KPS5wM7CxZ7kzaeu0TEd3AAWBCyrFV5St27ETgwM++cTpMw7olw/45aUK/tx+38j8P+uqTZiySWiS1SWrbs2dPipIqx1fsmNnx4swY/lmHNKHfCTSULNcDXX31SaZ3TgP2pxxLRCyLiEJEFOrq6tJXXwG9XcljZlYNuzVx2D8jTeivAaZKapJUC8wDWsv6tALzk9dXA49H8WRBKzAvubqnCZgK/LoypVdIjyt5iidxzY43/Zx6sww4FLXsvPCWYf+cfm+4FhHdkhYCq4Aa4P6I2CJpMdAWEa3Ad4GHJHVQPMKfl4zdImkFsBXoBm56vSt3qqb5Wmi+FgFtrffSsG4JZ8ZeDugNFM+u/z+/9uuqvd6tifxmwrGrd/xzmcXXuzWRnRcN7eqdtPq9emekjfTVO2ZmWVDJq3fMzCwjHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjDn0zsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY4cdzdck7QH+O0Q3mIikKfnH3p7sy1v2wv52+ZKbe+bI6LfB5Icd6E/VJLa0txpLiu8vdmWt+2F/G3zSG+vp3fMzHLEoW9mliNZDP1l1S5ghHl7sy1v2wv52+YR3d7MzembmVnfsnikb2ZmfchM6EuaLaldUoekRdWup9IkNUj6paTnJG2R9Nmk/QxJv5D0fPJ9fLVrrSRJNZLWS/pJstwk6Zlkex+RVFvtGitJ0umSHpX078m+/sss72NJn0t+njdLeljS2KztY0n3S9otaXNJW6/7VEXfSnJso6QLK11PJkJfUg2wFLgSmAZcJ2ladauquG7gf0bEucDFwE3JNi4CVkfEVGB1spwlnwWeK1n+O+CuZHv/CHyiKlUNn28CP4+ItwBvo7jtmdzHkiYDNwOFiDgPqAHmkb19/H+A2WVtfe3TK4GpyVcL8O1KF5OJ0AdmAB0RsS0iDgPLgblVrqmiIuJ3EbEuef2fFMNgMsXtfDDp9iBwVXUqrDxJ9cD7gPuSZQHvBh5NumRte08FLgO+CxARhyPiRTK8j4HRwDhJo4GTgd+RsX0cEb8C9pc197VP5wL/GEVPA6dL+q+VrCcroT8Z2Fmy3Jm0ZZKkKcB04Bngv0TE76D4iwE4s3qVVdw3gC8AryTLE4AXI6I7Wc7afj4L2AM8kExp3SfpDWR0H0fELuDvgR0Uw/4AsJZs7+Nj+tqnw55lWQl99dKWycuSJJ0C/Aj4m4j4U7XrGS6S3g/sjoi1pc29dM3Sfh4NXAh8OyKmA38mI1M5vUnmsecCTcAk4A0UpzfKZWkf92fYf8azEvqdQEPJcj3QVaVaho2kMRQD//sR8eOk+Q/H/vxLvu+uVn0VdgkwR9J2itN176Z45H96MhUA2dvPnUBnRDyTLD9K8ZdAVvfxFcBvImJPRBwBfgy8k2zv42P62qfDnmVZCf01wNTkrH8txZNBrVWuqaKS+ezvAs9FxJ0lq1qB+cnr+cA/jXRtwyEibo2I+oiYQnF/Ph4RHwZ+CVyddMvM9gJExO+BnZLOSZpmAlvJ6D6mOK1zsaSTk5/vY9ub2X1coq992gp8NLmK52LgwLFpoIqJiEx8AX8F/AfwAnBbtesZhu27lOKfeRuBDcnXX1Gc514NPJ98P6PatQ7Dtl8O/CR5fRbwa6AD+CFwUrXrq/C2XgC0Jft5JTA+y/sY+N/AvwObgYeAk7K2j4GHKZ6zOELxSP4Tfe1TitM7S5Mc20TxyqaK1uP/yDUzy5GsTO+YmVkKDn0zsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY449M3McuT/A1UioZaAkOmRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c8846cc908>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.misc import factorial\n",
    "n=100\n",
    "p=0.01\n",
    "x=np.arange(1,101)\n",
    "pois=1/factorial(x)/math.e\n",
    "bino=nCr(n,x)*np.power(p,x)*np.power(1-p,n-x)\n",
    "plt.scatter(x,pois)\n",
    "plt.scatter(x,bino)\n",
    "plt.legend([\"pois\",\"bino\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# problem10(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 1, 2, 0, 1, 0, 3, 1, 1, 3, 0, 0, 0, 0, 3, 0, 4, 0, 1, 3,\n",
       "       2, 0, 1, 2, 2, 1, 0, 0, 1, 1, 2, 1, 0, 0, 0, 0, 0, 1, 1, 2, 3, 0,\n",
       "       1, 2, 1, 0, 2, 0, 2, 0, 1, 0, 1, 2, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1,\n",
       "       1, 0, 0, 2, 0, 2, 1, 1, 0, 0, 3, 0, 1, 0, 2, 0, 0, 1, 2, 1, 3, 0,\n",
       "       0, 2, 2, 2, 1, 2, 2, 1, 4, 1, 1, 2, 1, 0, 1, 0, 2, 0, 1, 2, 1, 0,\n",
       "       1, 1, 1, 1, 2, 0, 0, 2, 3, 0, 1, 2, 2, 1, 0, 0, 0, 0, 2, 3, 1, 2,\n",
       "       1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 2, 0, 1, 1, 0, 0, 2, 1, 0,\n",
       "       0, 0, 1, 2, 3, 2, 0, 1, 0, 0, 1, 1, 2, 0, 0, 0, 2, 0, 1, 1, 0, 1,\n",
       "       1, 1, 0, 2, 1, 1, 0, 1, 1, 1, 1, 1, 3, 2, 1, 0, 0, 2, 2, 1, 0, 0,\n",
       "       4, 0, 0, 0, 1, 1, 0, 0, 5, 0, 0, 2, 2, 1, 2, 2, 0, 1, 0, 1, 2, 2,\n",
       "       1, 4, 2, 1, 0, 3, 0, 1, 2, 0, 1, 3, 0, 2, 2, 0, 1, 2, 2, 2, 0, 2,\n",
       "       0, 0, 0, 1, 1, 1, 0, 0, 2, 1, 1, 0, 1, 1, 2, 1, 0, 0, 0, 1, 1, 1,\n",
       "       2, 1, 3, 1, 0, 1, 2, 0, 1, 2, 1, 1, 1, 0, 1, 0, 4, 0, 1, 1, 2, 2,\n",
       "       1, 0, 0, 1, 1, 2, 1, 1, 1, 1, 1, 0, 2, 1, 1, 2, 3, 1, 1, 0, 1, 0,\n",
       "       0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 0, 2, 1, 2, 3, 1, 0,\n",
       "       4, 1, 2, 2, 0, 1, 2, 0, 0, 1, 1, 1, 1, 2, 0, 0, 1, 0, 1, 1, 2, 3,\n",
       "       0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 1, 1, 0, 0, 2, 0,\n",
       "       1, 0, 2, 0, 3, 0, 1, 1, 3, 1, 0, 1, 1, 2, 3, 2, 2, 2, 0, 1, 2, 3,\n",
       "       1, 0, 1, 1, 0, 1, 1, 2, 0, 1, 2, 0, 1, 1, 0, 0, 1, 2, 1, 0, 1, 3,\n",
       "       2, 2, 3, 1, 0, 0, 3, 1, 0, 1, 1, 1, 1, 2, 2, 0, 1, 0, 3, 3, 2, 0,\n",
       "       1, 2, 0, 2, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 3, 1, 0, 2, 1, 2, 0, 1,\n",
       "       2, 1, 0, 0, 0, 4, 0, 1, 1, 0, 0, 0, 0, 1, 0, 3, 1, 2, 2, 2, 0, 0,\n",
       "       1, 2, 3, 0, 3, 2, 2, 1, 2, 1, 2, 1, 0, 0, 2, 0, 1, 0, 0, 0, 1, 1,\n",
       "       1, 1, 2, 0, 0, 1, 2, 1, 1, 1, 1, 0, 1, 3, 2, 2, 1, 1, 0, 2, 1, 1,\n",
       "       1, 0, 0, 0, 0, 0, 2, 0, 1, 1, 0, 4, 3, 0, 1, 3, 0, 1, 3, 1, 0, 1,\n",
       "       1, 2, 0, 2, 1, 0, 3, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 2, 0, 2, 2,\n",
       "       0, 0, 0, 2, 2, 0, 0, 1, 0, 2, 2, 2, 0, 0, 1, 5, 0, 3, 0, 1, 0, 2,\n",
       "       0, 1, 1, 1, 1, 1, 2, 0, 1, 1, 0, 0, 2, 0, 0, 1, 1, 2, 0, 0, 1, 2,\n",
       "       0, 0, 1, 1, 0, 2, 1, 1, 1, 1, 2, 0, 2, 2, 1, 1, 0, 0, 1, 0, 0, 0,\n",
       "       1, 0, 2, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 3, 0, 2, 1, 4, 0, 3, 1,\n",
       "       1, 1, 2, 2, 3, 0, 0, 0, 0, 1, 1, 0, 2, 0, 1, 2, 0, 2, 1, 0, 1, 1,\n",
       "       0, 1, 0, 2, 0, 1, 1, 0, 0, 2, 0, 2, 0, 1, 0, 2, 3, 0, 2, 1, 1, 2,\n",
       "       1, 0, 1, 2, 2, 0, 1, 0, 2, 0, 1, 1, 1, 1, 0, 3, 0, 1, 1, 1, 2, 0,\n",
       "       3, 1, 0, 0, 1, 1, 2, 1, 1, 2, 1, 1, 0, 0, 5, 1, 0, 0, 1, 2, 3, 0,\n",
       "       1, 1, 1, 0, 0, 1, 1, 0, 0, 2, 1, 1, 0, 1, 1, 1, 2, 3, 1, 2, 2, 1,\n",
       "       1, 0, 1, 1, 1, 1, 0, 0, 2, 2, 1, 2, 3, 0, 1, 1, 4, 2, 0, 1, 0, 1,\n",
       "       0, 0, 1, 4, 2, 1, 2, 3, 2, 0, 2, 0, 1, 3, 0, 1, 2, 3, 1, 1, 0, 0,\n",
       "       1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 2, 1, 2, 1, 0, 1, 1, 2, 1, 3, 0, 0,\n",
       "       2, 0, 3, 1, 2, 0, 0, 1, 3, 1, 0, 1, 1, 2, 0, 0, 3, 1, 1, 0, 0, 2,\n",
       "       1, 1, 1, 0, 0, 0, 0, 1, 2, 5, 0, 1, 2, 0, 0, 0, 0, 0, 1, 1, 1, 1,\n",
       "       1, 0, 0, 1, 0, 0, 1, 2, 2, 0, 0, 1, 1, 2, 0, 1, 1, 0, 0, 1, 1, 3,\n",
       "       0, 3, 0, 1, 1, 1, 0, 2, 1, 0, 1, 0, 0, 1, 3, 1, 1, 0, 1, 2, 3, 1,\n",
       "       0, 0, 2, 1, 1, 0, 0, 1, 2, 5, 1, 1, 0, 2, 1, 1, 3, 1, 1, 2, 0, 1,\n",
       "       0, 0, 2, 0, 0, 0, 1, 0, 1, 1, 2, 1, 2, 0, 0, 0, 2, 0, 1, 2, 2, 0,\n",
       "       1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 2, 1, 0, 4, 0,\n",
       "       2, 2, 2, 3, 0, 2, 1, 3, 1, 1])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mock_possion_sample=np.random.binomial(100,0.01,1000)\n",
    "mock_possion_sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mathematical result: \n",
    "\n",
    "$ P(X>\\mu+a\\sigma) = P(\\frac{X-\\mu}{\\sigma}>a)=1-\\Psi(a)=\\Psi(-a)$\n",
    "\n",
    "$\\Psi(-1)=0.1587$\n",
    "\n",
    "$\\Psi(-2)=0.0228$\n",
    "\n",
    "$\\Psi(-2)=0.0013$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "when a= 1  p(x>miu+a*sigma)= 0.162800\n",
      "when a= 2  p(x>miu+a*sigma)= 0.022800\n",
      "when a= 3  p(x>miu+a*sigma)= 0.001100\n"
     ]
    }
   ],
   "source": [
    "miu=0\n",
    "sigma=1\n",
    "n=10000\n",
    "sample=np.random.normal(loc=miu,scale=sigma,size=n)\n",
    "for i in range(1,4):\n",
    "    p=sum(sample>miu+i*sigma)/n\n",
    "    print(\"when a= %d  p(x>miu+a*sigma)= %f\" %(i,p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the simulation for expected value of geometric distribution when p=0.3 is 3.323000, its variance is 7.720671\n"
     ]
    }
   ],
   "source": [
    "x=[]\n",
    "p=0.3\n",
    "for i in range(1000):\n",
    "    count=0\n",
    "    i=1\n",
    "    while i >p:\n",
    "        count=count+1\n",
    "        i=np.random.uniform()\n",
    "    x.append(count)\n",
    "print(\"the simulation for expected value of geometric distribution when p=0.3 is %f, its variance is %f\" %(np.mean(x),np.var(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The theoretical expect of geometic distribution when p=0.3 is $\\frac{1}{p}=3.33$\n",
    "\n",
    "The theoretical variance of geometic distribution when p=0.3 is $\\frac{1-p}{p^2}=7.78$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the simulation for expected value of exponential distribution when lamda= 2 is 0.468689, its variance is 0.207164\n"
     ]
    }
   ],
   "source": [
    "U=np.random.uniform(0,1,1000)\n",
    "lamda=2\n",
    "x=-1/lamda*np.log(U)\n",
    "print(\"the simulation for expected value of exponential distribution when lamda= 2 is %f, its variance is %f\" %(np.mean(x),np.var(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The theoretical expect of exponential distribution when p=0.3 is $\\frac{1}{\\lambda}=0.5$\n",
    "\n",
    "The theoretical variance of exponential distribution when p=0.3 is $\\frac{1}{\\lambda^2}=0.25$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}