Statistical_simulations_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistial Simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Libraries\n",
    "\n",
    "Python, like other programming languages, has an abundance of additional modules or libraries that augument the base framework and functionality of the language.\n",
    "\n",
    "Think of a library as a collection of functions that can be accessed to complete certain programming tasks without having to write your own algorithm.\n",
    "\n",
    "For this course, we will focus primarily on the following libraries:\n",
    "\n",
    "* **Numpy** is a library for working with arrays of data.\n",
    "\n",
    "* **Pandas** provides high-performance, easy-to-use data structures and data analysis tools.\n",
    "\n",
    "* **Scipy** is a library of techniques for numerical and scientific computing.\n",
    "\n",
    "* **Matplotlib** is a library for making graphs.\n",
    "\n",
    "* **Seaborn** is a higher-level interface to Matplotlib that can be used to simplify many graphing tasks.\n",
    "\n",
    "* **Statsmodels** is a library that implements many statistical techniques.\n",
    "\n",
    "# Documentation\n",
    "\n",
    "Reliable and accesible documentation is an absolute necessity when it comes to knowledge transfer of programming languages.  Luckily, python provides a significant amount of detailed documentation that explains the ins and outs of the language syntax, libraries, and more.  \n",
    "\n",
    "Understanding how to read documentation is crucial for any programmer as it will serve as a fantastic resource when learning the intricacies of python.\n",
    "\n",
    "Here is the link to the documentation of the python standard library: [Python Standard Library](https://docs.python.org/3/library/index.html#library-index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sn\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import statsmodels as stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.ads bid 如果有两个bidder X AND Y， the bid amount is uniform[0,1] distribution. what is the expected revenue?\n",
    "what if three bidders?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# two bidders\n",
    "x = np.random.uniform(low=0,high=1,size=10000)\n",
    "y = np.random.uniform(0,1,size=10000)\n",
    "z = np.random.uniform(0,1,size=10000)\n",
    "xyz = np.stack((x,y,z),axis=0)\n",
    "xyz_sort = np.sort(xyz,axis=0)\n",
    "exp1 = np.mean(np.maximum(x,y))## maximum take two ndarrys where max take one ndarray\n",
    "exp2 = np.mean(np.minimum(x,y))\n",
    "exp3 = np.mean(np.minimum(x,y,z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 非常简单的for loop coding 题。用R或者python写都可以。计算 sqrt(1) + sqrt(2) + ... + sqrt(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "671.4629471031477"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math as mt\n",
    "sum([mt.sqrt(i) for i in range(1,101)])\n",
    "sum([np.sqrt(i) for i in range(1,101)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. generate poisson distribution using unifrom\n",
    "https://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def poisson(lam):\n",
    "    k, p = 0, 1\n",
    "    while p > mt.exp(-lam):\n",
    "        k += 1\n",
    "        p = p*np.random.uniform()\n",
    "    return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poisson(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bootstrapped Mean Length = 0.4909181528776206, 95% CI = [0.47222953 0.50840311]\n"
     ]
    }
   ],
   "source": [
    "# bootstrapping\n",
    "# Draw some random sample with replacement and append mean to mean_lengths.\n",
    "mean_lengths, sims = [], 1000\n",
    "wrench_lengths = np.random.uniform(size=1000)\n",
    "for i in range(sims):\n",
    "    temp_sample = np.random.choice(wrench_lengths, replace=True, size=len(wrench_lengths))\n",
    "    sample_mean = np.mean(temp_sample)\n",
    "    mean_lengths.append(sample_mean)\n",
    "    \n",
    "# Calculate bootstrapped mean and 95% confidence interval.\n",
    "boot_mean = np.mean(mean_lengths)\n",
    "boot_95_ci = np.percentile(mean_lengths, [2.5, 97.5])\n",
    "print(\"Bootstrapped Mean Length = {}, 95% CI = {}\".format(boot_mean, boot_95_ci))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Non-standard estimators\n",
    "In the last exercise, you ran a simple bootstrap that we will now modify for more complicated estimators.\n",
    "\n",
    "Suppose you are studying the health of students. You are given the height and weight of 1000 students and are interested in the median height as well as the correlation between height and weight and the associated 95% CI for these quantities. Let's use bootstrapping.\n",
    "\n",
    "Examine the pandas DataFrame df with the heights and weights of 1000 students. Using this, calculate the 95% CI for both the median height as well as the correlation between height and weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['ID', 'Age', 'Gender', 'GenderGroup', 'Glasses', 'GlassesGroup',\n",
       "       'Height', 'Wingspan', 'CWDistance', 'Complete', 'CompleteGroup',\n",
       "       'Score'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Height Median CI = [65. 70.] \n",
      "Height Weight Correlation CI = [-0.55614018  1.        ]\n"
     ]
    }
   ],
   "source": [
    "# Sample with replacement and calculate quantities of interest\n",
    "# Store the url string that hosts our .csv file\n",
    "url = \"Cartwheeldata.csv\"\n",
    "\n",
    "# Read the .csv file and store it as a pandas Data Frame\n",
    "df = pd.read_csv(url)\n",
    "\n",
    "sims, data_size, height_medians, hw_corr = 1000, df.shape[0], [], []\n",
    "for i in range(sims):\n",
    "    tmp_df = df.sample(n=data_size,replace=True)\n",
    "    height_medians.append(tmp_df.median().Height)\n",
    "    hw_corr.append(tmp_df.corr().Height)\n",
    "\n",
    "# Calculate confidence intervals\n",
    "height_median_ci = np.percentile(height_medians,[2.5,97.5])\n",
    "height_weight_corr_ci = np.percentile(hw_corr,[2.5,97.5])\n",
    "print(\"Height Median CI = {} \\nHeight Weight Correlation CI = {}\".format( height_median_ci, height_weight_corr_ci))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bootstrapping regression\n",
    "Now let's see how bootstrapping works with regression. Bootstrapping helps estimate the uncertainty of non-standard estimators. Consider the R2 statistic associated with a regression. When you run a simple least squares regression, you get a value for R2. But let's see how can we get a 95% CI for R2.\n",
    "\n",
    "Examine the DataFrame df with a dependent variable y and two independent variables X1 and X2 using df.head(). We've already fit this regression with statsmodels (sm) using:\n",
    "\n",
    "reg_fit = sm.OLS(df['y'], df.iloc[:,1:]).fit()\n",
    "Examine the result using reg_fit.summary() to find that R2=0.3504. Use bootstrapping to calculate the 95% CI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1.])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "rsquared_boot, coefs_boot, sims = [], [], 1000\n",
    "reg_fit = sm.OLS(df['Height'], df.iloc[:,6]).fit()\n",
    "\n",
    "# Run 1K iterations\n",
    "for i in range(sims):\n",
    "    # First create a bootstrap sample with replacement with n=df.shape[0]\n",
    "    bootstrap = df.sample(n=df.shape[0],replace=True)\n",
    "    # Fit the regression and append the r square to rsquared_boot\n",
    "    rsquared_boot.append(sm.OLS(bootstrap['Height'],bootstrap.iloc[:,6]).fit().rsquared)\n",
    "\n",
    "# Calculate 95% CI on rsquared_boot\n",
    "r_sq_95_ci = np.percentile(rsquared_boot,[2.5,97.5])\n",
    "r_sq_95_ci"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next few exercises, we will run a significance test using permutation testing. \n",
    "Hypothesis testing - Difference of means\n",
    "We want to test the hypothesis that there is a difference in the average donations received from A and B. Previously, you learned how to generate one permutation of the data. Now, we will generate a null distribution of the difference in means and then calculate the p-value.\n",
    "\n",
    "For the null distribution, we first generate multiple permuted datasets and store the difference in means for each case. We then calculate the test statistic as the difference in means with the original dataset. Finally, we calculate the p-value as twice the fraction of cases where the difference is greater than or equal to the absolute value of the test statistic (2-sided hypothesis). A p-value of less than say 0.05 could then determine statistical significance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value = 0.78\n"
     ]
    }
   ],
   "source": [
    "# Generate permutations equal to the number of repetitions\n",
    "data = np.random.uniform(size = 2000)\n",
    "donations_A = range(100)\n",
    "donations_B = range(100)\n",
    "reps = 100\n",
    "perm = np.array([np.random.permutation(len(donations_A) + len(donations_B)) for i in range(reps)])\n",
    "permuted_A_datasets = data[perm[:, :len(donations_A)]]\n",
    "permuted_B_datasets = data[perm[:, len(donations_A):]]\n",
    "\n",
    "# Calculate the difference in means for each of the datasets\n",
    "samples = np.mean(permuted_A_datasets,axis=1) - np.mean(permuted_B_datasets,axis=1)\n",
    "\n",
    "# Calculate the test statistic and p-value\n",
    "test_stat = np.mean(donations_A) - np.mean(donations_B)\n",
    "p_val = 2*np.sum(samples >= np.abs(test_stat))/reps\n",
    "print(\"p-value = {}\".format(p_val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(permuted_A_datasets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hypothesis testing - Non-standard statistics\n",
    "In the previous two exercises, we ran a permutation test for the difference in mean values. Now let's look at non-standard statistics.\n",
    "\n",
    "Suppose that you're interested in understanding the distribution of the donations received from websites A and B. For this, you want to see if there's a statistically significant difference in the median and the 80th percentile of the donations. Permutation testing gives you a wonderfully flexible framework for attacking such problems.\n",
    "\n",
    "Let's go through running a test to see if there's a difference in the median and the 80th percentile of the distribution of donations. As before, you're given the donations from the websites A and B in the variables donations_A and donations_B respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "80th Percentile: test statistic = 0.0, p-value = 0.88\n",
      "Median: test statistic = 0.0, p-value = 0.88\n"
     ]
    }
   ],
   "source": [
    "# Calculate the difference in 80th percentile and median for each of the permuted datasets (A and B)\n",
    "samples_percentile = np.percentile(permuted_A_datasets, 80, axis=1) - np.percentile(permuted_B_datasets, 80, axis=1)\n",
    "samples_median = np.median(permuted_A_datasets, axis=1) - np.median(permuted_B_datasets, axis=1)\n",
    "\n",
    "# Calculate the test statistic from the original dataset and corresponding p-values\n",
    "test_stat_percentile = np.percentile(donations_A, 80) - np.percentile(donations_B, 80)\n",
    "test_stat_median = np.median(donations_A) - np.median(donations_B)\n",
    "p_val_percentile = 2*np.sum(samples_percentile >= np.abs(test_stat_percentile))/reps\n",
    "p_val_median = 2*np.sum(samples_median >= np.abs(test_stat_median))/reps\n",
    "\n",
    "print(\"80th Percentile: test statistic = {}, p-value = {}\".format(test_stat_percentile, p_val_percentile))\n",
    "print(\"Median: test statistic = {}, p-value = {}\".format(test_stat_median, p_val_median))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Power Analysis - Part II\n",
    "Previously, we simulated one instance of the experiment & generated a p-value. We will now use this framework to calculate statistical power. Power of an experiment is the experiment's ability to detect a difference between treatment & control if the difference really exists. It's good statistical hygiene to strive for 80% power.\n",
    "\n",
    "For our website, we want to know how many people need to visit each variant, such that we can detect a 10% increase in time spent with 80% power. For this, we start with a small sample (50), simulate multiple instances of this experiment & check power. If 80% power is reached, we stop. If not, we increase the sample size & try again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# statsmodels\n",
    "https://www.statsmodels.org/stable/_modules/statsmodels/stats/power.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For 80% power, sample size required = 1510\n"
     ]
    }
   ],
   "source": [
    "import scipy.stats as st\n",
    "sample_size = 50\n",
    "control_mean = 1\n",
    "control_sd = 1\n",
    "effect_size = 0.1\n",
    "sim = 100\n",
    "# Keep incrementing sample size by 10 till we reach required power\n",
    "while 1:\n",
    "    control_time_spent = np.random.normal(loc=control_mean, scale=control_sd, size=(sample_size,sims))\n",
    "    treatment_time_spent = np.random.normal(loc=control_mean*(1+effect_size), scale=control_sd, size=(sample_size,sims))\n",
    "    t, p = st.ttest_ind(treatment_time_spent, control_time_spent)\n",
    "    \n",
    "    # Power is the fraction of times in the simulation when the p-value was less than 0.05\n",
    "    power = (p < 0.05).sum()/sims\n",
    "    if power > 0.8: \n",
    "        break;\n",
    "    else: \n",
    "        sample_size += 10\n",
    "print(\"For 80% power, sample size required = {}\".format(sample_size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a simple exercise introducing the concept of Monte Carlo Integration.\n",
    "\n",
    "Here we will evaluate a simple integral ∫10xexdx. We know that the exact answer is 1, but simulation will give us an approximate solution, so we can expect an answer close to 1. As we saw in the video, it's a simple process. For a function of a single variable f(x):\n",
    "\n",
    "Get the limits of the x-axis (xmin,xmax) and y-axis (max(f(x)),min(min(f(x)),0)).\n",
    "Generate a number of uniformly distributed point in this box.\n",
    "Multiply the area of the box ((max(f(x)−min(f(x))×(xmax−xmin)) by the fraction of points that lie below f(x).\n",
    "Upon completion, you will have a framework for handling definite integrals using Monte Carlo Integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated answer = 0.9788947624753406, Actual Answer = 1\n"
     ]
    }
   ],
   "source": [
    "# Define the sim_integrate function\n",
    "def sim_integrate(func, xmin, xmax, sims):\n",
    "    x = np.random.uniform(xmin, xmax, sims)\n",
    "    y = np.random.uniform(min(min(func(x)),0), max(func(x)), sims)\n",
    "    area = (max(y) - min(y))*(xmax-xmin)\n",
    "    result = area * sum( (y-np.array(func(x)))<=0 )/sims\n",
    "    return result\n",
    "    \n",
    "# Call the sim_integrate function and print results\n",
    "result = sim_integrate(func = lambda x: x*np.exp(x), xmin = 0, xmax = 1, sims = 500)\n",
    "print(\"Simulated answer = {}, Actual Answer = 1\".format(result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the value of pi\n",
    "Now we work through a classic example - estimating the value of π.\n",
    "\n",
    "Imagine a square of side 2 with the origin (0,0) as its center and the four corners having coordinates (1,1),(1,−1),(−1,1),(−1,−1). The area of this square is 2×2=4. Now imagine a circle of radius 1 with its center at the origin fitting perfectly inside this square. The area of the circle will be π×radius2=π.\n",
    "\n",
    "To estimate π, we randomly sample multiple points in this square & get the fraction of points inside the circle (x2+y2<=1). The area of the circle then is 4 times this fraction, which gives us our estimate of π.\n",
    "\n",
    "After this exercise, you'll have a grasp of how to use simulation for computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated value of pi = 3.1284\n"
     ]
    }
   ],
   "source": [
    "# Initialize sims and circle_points\n",
    "sims, circle_points = 10000, 0\n",
    "\n",
    "for i in range(sims):\n",
    "    # Generate the two coordinates of a point\n",
    "    point = np.random.uniform(-1,1,size=2)\n",
    "    # if the point lies within the unit circle, increment counter\n",
    "    within_circle = point[0]**2 + point[1]**2 <= 1\n",
    "    if within_circle == True:\n",
    "        circle_points +=1\n",
    "        \n",
    "# Estimate pi as 4 times the avg number of points in the circle.\n",
    "pi_sim = 4*circle_points/sims\n",
    "print(\"Simulated value of pi = {}\".format(pi_sim))"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}