MLE and KS Test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assignment 4\n",
    "\n",
    "By Siddharth Sundararajan and Nitin Reddy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing package into 'C:/Users/nitin/OneDrive/Documents/R/win-library/3.5'\n",
      "(as 'lib' is unspecified)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'matlib' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\nitin\\AppData\\Local\\Temp\\RtmpyIYYzl\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "# Installing the required packages\n",
    "install.packages(\"matlib\")\n",
    "library(matlib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Writing the required symbols\n",
    "mu=intToUtf8(0x03BC)\n",
    "sigma=intToUtf8(0x03C3)\n",
    "varcovar = intToUtf8(0x03A3)\n",
    "alpha = intToUtf8(0x03B1)\n",
    "beta = intToUtf8(0x03B2)\n",
    "lambda = intToUtf8(0x03BB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MLE for all the functions\n",
    "\n",
    "# Bernoulli function\n",
    "bernoulli.mle <- function(x) {\n",
    "  p <- mean(x)\n",
    "  cat(\"BERNOULLI DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for bernoulli distribution is : p = \", p, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(p)\n",
    "}\n",
    "\n",
    "# Binomial function\n",
    "binomial.mle <- function(x) {\n",
    "  p <- mean(x)\n",
    "  n <- length(x)\n",
    "  cat(\"BINOMIAL DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for binomial distribution is : n = \", n,\n",
    "      \", p = \", p, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(c(n,p))\n",
    "}\n",
    "\n",
    "# Geometric function\n",
    "geometric.mle <- function(x) {\n",
    "  p <- 1/(mean(x))\n",
    "  cat(\"GEOMETRIC DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for geometric distribution is : p = \", p, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(p)\n",
    "}\n",
    "\n",
    "# Poisson function\n",
    "poisson.mle <- function(x) {\n",
    "  lambda.val <- mean(x)\n",
    "  cat(\"POISSON DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for poisson distribution is : \", lambda,\n",
    "      \" = \", lambda.val, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(lambda.val)\n",
    "}\n",
    "\n",
    "# Uniform function\n",
    "uniform.mle <- function(x) {\n",
    "  a <- min(x)\n",
    "  b <- max(x)\n",
    "  cat(\"UNIFORM DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for uniform distribution are : a = \", a,\n",
    "      \", b = \", b, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(c(a,b))\n",
    "}\n",
    "\n",
    "# Normal function\n",
    "normal.mle <- function(x) {\n",
    "  mu.val <- mean(x)\n",
    "  sigma.val <- sum(((x - mean(x))^2))/length(x)\n",
    "  cat(\"NORMAL DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for normal distribution are : \", mu, \" = \",\n",
    "      mu.val, \", \", sigma, \" = \", sigma.val, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(c(mu.val,sigma.val))\n",
    "}\n",
    "\n",
    "# Exponential function\n",
    "exponential.mle <- function(x) {\n",
    "  beta.val <- mean(x)\n",
    "  cat(\"EXPONENTIAL DISTRIBUTION \\n\")\n",
    "  cat(\"The population estimate for exponential distribution is : \", beta,\n",
    "      \" = \", beta.val, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(beta.val)\n",
    "}\n",
    "\n",
    "# Gamma function\n",
    "gamma.mle <- function(x) {\n",
    "  cat(\"GAMMA DISTRIBUTION \\n\")\n",
    "  mean.val <- mean(x)\n",
    "  s <- log(mean.val) - ((sum(log(x)))/length(x))\n",
    "  alpha.val <- (3 - s + sqrt( ((s-3)^2) + (24 * s) ))/(12 * s)\n",
    "  beta.val <- mean.val/alpha.val\n",
    "  cat(\"The population estimate for gamma distribution is : \", alpha, \" = \",\n",
    "      alpha.val, \", \", beta, \" = \", beta.val, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(c(alpha.val,beta.val))\n",
    "}\n",
    "\n",
    "# Beta function\n",
    "beta.mle <- function(x) {\n",
    "  cat(\"BETA DISTRIBUTION \\n\")\n",
    "  iter = 100\n",
    "  # Initial values for alpha and beta are from the MOM estimator\n",
    "  x.mean <- mean(x)\n",
    "  x.var <- (sum(x * x))/length(x)\n",
    "  alpha.val <- ((x.mean ^ 2) - (x.mean * x.var))/(x.var - (x.mean ^ 2))\n",
    "  beta.val <- (alpha.val * (1 - x.mean))/(x.mean)\n",
    "  # Numerical approximation\n",
    "  final.val <- c(alpha.val, beta.val)\n",
    "  for (i in 1:100) {\n",
    "    g1 <- digamma(alpha.val) - digamma(alpha.val + beta.val)\n",
    "      - (sum(log(x))/length(x))\n",
    "    g2 <- digamma(beta.val) - digamma(alpha.val + beta.val)\n",
    "      - (sum(log(1-x))/length(x))\n",
    "    g.val <- c(g1, g2)\n",
    "    # Matrix G terms\n",
    "    term.11 <- trigamma(alpha.val) - trigamma(alpha.val + beta.val)\n",
    "    term.12 <- -trigamma(alpha.val + beta.val)\n",
    "    term.22 <- trigamma(beta.val) - trigamma(alpha.val + beta.val)\n",
    "    G <- matrix(c(term.11, term.12, term.12, term.22),\n",
    "                nrow = 2, ncol = 2, byrow = TRUE)\n",
    "    G.inv <- inv(G)\n",
    "    # Final values\n",
    "    final.val <- final.val - t(G.inv %*% g.val)\n",
    "    alpha.val <- final.val[1]\n",
    "    beta.val <- final.val[2]\n",
    "  }\n",
    "  cat(\"The population estimate for beta distribution is : \", alpha, \" = \",\n",
    "      alpha.val, \", \", beta, \" = \", beta.val, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(c(c(alpha.val,beta.val)))\n",
    "}\n",
    "\n",
    "# Chisquare function\n",
    "chisquare.mle <- function(x) {\n",
    "  cat(\"CHISQUARE DISTRIBUTION \\n\")\n",
    "  # Intitial values for v from MOM estimator\n",
    "  p.tilda <- mean(x)\n",
    "  # Some approximation using second derivative\n",
    "  n <- length(x)\n",
    "  p.del.numerator <- (-n/gamma(p.tilda/2) * digamma(p.tilda/2))\n",
    "    - (((n * log(2)) + sum(log(x)))/2)\n",
    "  p.del.denominator <- (-n * trigamma(p.tilda/2)/4)\n",
    "  p.del <- p.del.numerator/p.del.denominator\n",
    "  # Not sure here\n",
    "  p.val <- (p.tilda + p.del)/2\n",
    "  cat(\"The population estimate for student-t distribution is : \n",
    "        p(degrees of freedom) = \", p.val, \"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  return(p.val)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enter distribution name\n",
    "cat(\"1. Bernoulli Distribution \\n\")\n",
    "cat(\"2. Binomial Distribution \\n\")\n",
    "cat(\"3. Geometric Distribution \\n\")\n",
    "cat(\"4. Poisson Distribution \\n\")\n",
    "cat(\"5. Uniform Distribution \\n\")\n",
    "cat(\"6. Normal Distribution \\n\")\n",
    "cat(\"7. Exponential Distribution \\n\")\n",
    "cat(\"8. Gamma Distribution \\n\")\n",
    "cat(\"9. Beta Distribution \\n\")\n",
    "cat(\"10. Chi-square Distribution \\n\")\n",
    "in.val <- as.integer(readline(prompt = \"Please select the distribution\n",
    "                    (Enter the number) \\n\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taking sample input data from the user or \n",
    "# asking user to input parameters to a R function\n",
    "cat(\"\\n Please select one option : \\n\")\n",
    "cat(\" 1. Enter data \\n\")\n",
    "cat(\" 2. Create data using function in R \\n\")\n",
    "in.val.2 <- as.integer(readline(prompt = \"Enter the number \\n\"))\n",
    "\n",
    "if(in.val.2 == 1) {\n",
    "  my.input <- readline(prompt \n",
    "                       = \"Please enter the sample data separated with commas \\n\")\n",
    "  x <- as.double(strsplit(my.input, \",\")[[1]])\n",
    "} else {\n",
    "  obs.input <- as.integer(readline(prompt \n",
    "                    = \"Please enter the number of observations(n) : \\n\"))\n",
    "  if(in.val == 1) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of p : \\n\"))\n",
    "    x <- rbinom(obs.input, 1, value)\n",
    "  }\n",
    "  else if(in.val == 2) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of p : \\n\"))\n",
    "    value2 <- as.double(readline(prompt = \"Please enter the value of n : \\n\"))\n",
    "    x <- rbinom(obs.input,n = value2,prob = value)\n",
    "  }\n",
    "  else if(in.val == 3) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of p : \\n\"))\n",
    "    x <- rgeom(obs.input, value)\n",
    "  }\n",
    "  else if(in.val == 4) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of lambda : \\n\"))\n",
    "    x <- rpois(obs.input, value)\n",
    "  }\n",
    "  else if(in.val == 5) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of a : \\n\"))\n",
    "    value2 <- as.double(readline(prompt = \"Please enter the value of b : \\n\"))\n",
    "    x <- runif(obs.input, value, value2)\n",
    "  }\n",
    "  else if(in.val == 6) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of mu : \\n\"))\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of sigma^2 : \\n\"))\n",
    "    x <- rnorm(obs.input, value, value2)\n",
    "  }\n",
    "  else if(in.val == 7) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of beta : \\n\"))\n",
    "    x <- rexp(obs.input, value)\n",
    "  }\n",
    "  else if(in.val == 8) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of alpha : \\n\"))\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of beta : \\n\"))\n",
    "    x <- rgamma(obs.input,shape = value, scale =value2)\n",
    "  }\n",
    "  else if(in.val == 9) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of alpha : \\n\"))\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of beta : \\n\"))\n",
    "    x <- rbeta(obs.input,shape1 = value, shape2 = value2)\n",
    "  }\n",
    "  else if(in.val == 10) {\n",
    "    value <- as.double(readline(prompt = \"Please enter the value of p : \\n\"))\n",
    "    x <- rchisq(obs.input, value)\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CHISQUARE DISTRIBUTION \n",
      "The population estimate for student-t distribution is : p(degrees of freedom) =  8.915059 \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "8.9150594891422"
      ],
      "text/latex": [
       "8.9150594891422"
      ],
      "text/markdown": [
       "8.9150594891422"
      ],
      "text/plain": [
       "[1] 8.915059"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# MLE based on input given by the user\n",
    "if (in.val == 1) {\n",
    "  bernoulli.mle(x)\n",
    "} else if (in.val == 2) {\n",
    "  binomial.mle(x)\n",
    "} else if (in.val == 3) {\n",
    "  geometric.mle(x)\n",
    "} else if (in.val == 4) {\n",
    "  poisson.mle(x)\n",
    "} else if (in.val == 5) {\n",
    "  uniform.mle(x)\n",
    "} else if (in.val == 6) {\n",
    "  normal.mle(x)\n",
    "} else if (in.val == 7) {\n",
    "  exponential.mle(x)\n",
    "} else if (in.val == 8) {\n",
    "  gamma.mle(x)\n",
    "} else if (in.val == 9) {\n",
    "  beta.mle(x)\n",
    "} else if (in.val == 10) {\n",
    "  chisquare.mle(x)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametric bootstrap with KS Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function for parametric bootstrap with ks test\n",
    "#************************************************#\n",
    "#nboot = number of iterations for bootstrap#\n",
    "#val = value of parameter 1#\n",
    "#val2 = value of parameter 2#\n",
    "#n = number of observations for the distribution#\n",
    "#in.val.n = distribution selected by the user#\n",
    "#************************************************#\n",
    "ks.test.bootstrap.mle <- function(nboot = 10000, val, val2 = 0, n = 10000, in.val.n) \n",
    "{\n",
    "  D.vec <- NULL\n",
    "  if (in.val.n == 1) {\n",
    "    func.mle <- bernoulli.mle\n",
    "    func <- rbinom\n",
    "    term <- \"rbinom\"\n",
    "  } else if (in.val.n == 2) {\n",
    "    func.mle <- binomial.mle\n",
    "    func <- rbinom\n",
    "    term <- \"rbinom\"\n",
    "  } else if (in.val.n == 3) {\n",
    "    func.mle <- geometric.mle\n",
    "    func <- rgeom\n",
    "    term <- \"rgeom\"\n",
    "  } else if (in.val.n == 4) {\n",
    "    func.mle <- poisson.mle\n",
    "    func <- rpois\n",
    "    term <- \"rpois\"\n",
    "  } else if (in.val.n == 5) {\n",
    "    func.mle <- uniform.mle\n",
    "    func <- runif\n",
    "    term <- \"runif\"\n",
    "  } else if (in.val.n == 6) {\n",
    "    func.mle <- normal.mle\n",
    "    func <- rnorm\n",
    "    term <- \"rnorm\"\n",
    "  } else if (in.val.n == 7) {\n",
    "    func.mle <- exponential.mle\n",
    "    func <- rexp\n",
    "    term <- \"rexp\"\n",
    "  } else if (in.val.n == 8) {\n",
    "    func.mle <- gamma.mle\n",
    "    func <- rgamma\n",
    "    term <- \"rgamma\"\n",
    "  } else if (in.val.n == 9) {\n",
    "    func.mle <- beta.mle\n",
    "    func <- rbeta\n",
    "    term <- \"rbeta\"\n",
    "  } else if (in.val.n == 10) {\n",
    "    func.mle <- chisquare.mle\n",
    "    func <- rchisq\n",
    "    term <- \"rchisq\"\n",
    "  }\n",
    "  for(i in 1:nboot) {\n",
    "    if (in.val.n == 3 | in.val.n == 4 | in.val.n == 7 | in.val.n == 10) {\n",
    "      vecb <- func(n, val)\n",
    "      est.val <- func.mle(vecb)\n",
    "      D.vec <- c(D.vec, unname(ks.test(vecb, term, est.val)$statistic))\n",
    "    } else {\n",
    "      if (in.val.n == 1) {val = 1}\n",
    "      vecb <- func(n, val, val2) \n",
    "      est.val <- func.mle(vecb)\n",
    "      if (in.val.n == 1) {\n",
    "        est.val[2] = est.val[1]\n",
    "        est.val[1] = 1}\n",
    "      D.vec <- c(D.vec, unname(ks.test(vecb, term,\n",
    "                                       est.val[1], est.val[2])$statistic))\n",
    "    }\n",
    "  }\n",
    "  list(D.values = D.vec)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. Bernoulli Distribution \n",
      "2. Binomial Distribution \n",
      "3. Geometric Distribution \n",
      "4. Poisson Distribution \n",
      "5. Uniform Distribution \n",
      "6. Normal Distribution \n",
      "7. Exponential Distribution \n",
      "8. Gamma Distribution \n",
      "9. Beta Distribution \n",
      "10. Chi-square Distribution \n",
      "Please select the distribution (Enter the number) \n",
      "10\n"
     ]
    }
   ],
   "source": [
    "# Enter distribution name for which goodness of fit test has to be done\n",
    "cat(\"1. Bernoulli Distribution \\n\")\n",
    "cat(\"2. Binomial Distribution \\n\")\n",
    "cat(\"3. Geometric Distribution \\n\")\n",
    "cat(\"4. Poisson Distribution \\n\")\n",
    "cat(\"5. Uniform Distribution \\n\")\n",
    "cat(\"6. Normal Distribution \\n\")\n",
    "cat(\"7. Exponential Distribution \\n\")\n",
    "cat(\"8. Gamma Distribution \\n\")\n",
    "cat(\"9. Beta Distribution \\n\")\n",
    "cat(\"10. Chi-square Distribution \\n\")\n",
    "in.val.bootstrap <- as.integer(readline(prompt =\n",
    "                \"Please select the distribution (Enter the number) \\n\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# D value or the statistic from the the function written earlier\n",
    "#************************************************#\n",
    "#nboot = number of iterations for bootstrap#\n",
    "#val = value of parameter 1#\n",
    "#val2 = value of parameter 2#\n",
    "#n = number of observations for the distribution#\n",
    "#in.val.n = distribution selected by the user#\n",
    "#************************************************#\n",
    "d.values <- as.numeric(unlist(ks.test.bootstrap.mle(nboot = 10000,\n",
    "                            val = 6,val2 = 0.6, n = 1000, in.val.bootstrap)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAgAElEQVR4nO3di3aizBZF4eIi3vH937a5qWAltgV7hcKa3zjjdLp7h0Jh/kal\nE3cDsJhbeweAb0BIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQY\nICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQY\nICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIyzjnph89/2Bs9yc7s8+c\n81f6eY98n87hJ9x3y3wU0jn7k7t536xNSCvhvlvmo5D+6BTNnbv88MeE9Be475bxQno/9Ef7\nMmt1QlqC+26Z3x6R6n3RfFQe+z9x97HTrv3q6zR8yrX5XXEYfeY1d1Xz0bFsPs6r6317h9zl\n59vtkLniPF1+sr3XEq67zOWHyR+f71/7NZ93/mGl0Vbuv9ZV5rJ+YnyrMEVIy/wS0jUb8ikm\nIRXDx2X3Gedh5PmZefcJ96nuVB8Gmsiqx589jLc37vWnzfeyRybZzV9pfIuGX++35PxyqzBF\nSMv8ElLzH/zmP9t1c6IeRmd4eT9v+5Kyx2/vn+naTzs0J2p9u1WTPpoTf9xgb7I9L6SXzfea\nrbaPXyfXPvR5K41v0fDrfSvZy63CFCEt48aGP+j/v/1iqG4eZB5/1p697lA3Xx+57mw+Nmdn\n+0v2/Mz2tG5fM7hOttT8aXPKu/zS/fJc+2V7L1/aDZs/ZZM/vvRf2/Vf2Xkr3W6vIfWt1X1C\n01uFMUJa5peQ2rP38VTofm7u7v8lr7qzuexP//aMv3/m6WXT/f+fJ788B1629xLSffOn6R83\n8dSvKbwJqWzHu9+Vr7cKY4S0zC8h7fs/GM665191J+Xt2v3B/emK99fNwLEq3E+n9ySKl+29\nhOReNj/Yt48tTbv7/6z0rOv5td30VmGMkJZ5nqXT87C6n35X76/uH72e6Y/fH/NRme9Dmnz0\nUUjX9qWCYvia7s1KXkju9VZhjJCW+S2kW33sXxIrJn/1eATJfnxE6n7bfqmX7w6XoEek7PUv\nfwupjegyfGX3bqXn16jjTx7fKowR0jK/htTq3uV5/ln53+dI3d/mw5//N6Tyo+dIx5eQmt/f\nP/G3ldo8z+7+HOn1C7n7rcIY98gyv4SUP56j3x8q6l9ftXM/PhTc/v+I9P5Vu8PLi4J3dfeQ\nUv+6Uta9Mn4ePq3dyXP3S/F6qzBGSMv8ElJzjhfX7tl5e6VCe1K2vz7e/+wvL/DfR+o2VHTD\n95et34T0ur2XYqabf/7l7vGl2Y8r7aZ79djK+fVWYYyQlvntS7v70/LujH2cuvczv79M5zSM\nvEQyXJHQvgN7fh/Sy/Ze/nLYTvkaUrvs8feVrv2fVcNvh53s05ncKowR0jK/hdQ/kyiGKwDK\nRzu7bPTy8aW91u7kRdL+cba7XO8XLPyw9bvJ9l7/cnIp3+gv3eN19p9W6vfq+HzVpGq+oitP\n9/VGtwojhLS6mqccX4CQVjN8hXUpphfQYZMIaTXPlwq8V5ixOYS0msc/SuA1sC9ASOup9+0r\nahmXrn0DQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYI\nCTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYI\nCTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYI\nCTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEJLHfWbt3URUOB88n90l3HEY\n43zwEBLCcT54CAnhOB88hIRwnA8eQkI4zgcPISEc54OHkBCO88FDSAjH+eAhJISbfz6c92X3\nDn9ZnQ33JwKEhHBzz4c6H10tU5ju0toICeHmng+Vy46X7qPrKXOV3Q6tj5AQbu75kLnL4+OL\ny2x2Jg6EhHBzz4fJ1c/fdSk0ISEcj0geQkK4Bc+RTtfuI54jAfPPh2L0ql1eW+7S2ggJ4Ra8\nj1R17yNl5Z73kZA8zgcPISEc54OHkBCOS4Q8hIRwXCLkISSE4xIhDyEhHG/IeggJ4bhEyENI\nCMcjkoeQEI5LhDyEhHBcIuQhJITjEiEPISEc54OHkBBOdD5s+QegEBLCLTwfDrlz5Um6xJ8j\nJIRb9j7S8IrD+xfttnbKERLCLQqpclV9u10rd1AssRZCQrhFIWWue927drliibUQEsItCun+\nOsL71xO2dsoREsItCml3D+ntJUJbO+UICeHmh1TuDyd3bD6sq/evNmztlCMkhJsf0uM9Iuey\nt5cIbe2UIySEm30+XC6HQ1l2LzlU7y+129opR0gI9wfnw9ZOOUJCOELyEBLCzT4f6p1zxXBx\nEC9/I3Wzv4tQ1n8vrn4jhITEzf8XsoempkPWfScuQkLq5n/Phu6Xa5ZfCQlY+l2E6qIgJGDu\n+ZC7+5tHeUFISN7c8+HgdsNHV1cQElI3+3yoHvWc/vOvybd2yhESws0/Hy7l/aPrjpCQOK5s\n8BASwhGSh5AQjpA8hIRwhOQhJIQjJA8hIRwheQgJ4QjJQ0gIR0geQkI4QvIQEsIRkoeQEI6Q\nPISEcITkISSEIyQPISEcIXkICeEIyUNICEdIHkJCOELyEBLCEZKHkBCOkDyEhHCE5CEkhCMk\nDyEhHCF5CAnhCMlDSAhHSB5CQjhC8hASwhGSh5AQjpA8hIRwhOQhJIQjJA8hIRwheQgJ4QjJ\nQ0gIR0geQkI4QvIQEsIRkufDkD4i3lVEg5A8ho9IW7vpmI2QPISEcITkISSEIyQPISEcIXkI\nCeEIyUNICEdIHkJCOELyEBLCEZKHkBCOkDyEhHCE5CEkhCMkDyEhHCF5CAnhCMlDSAhHSB5C\nQjhC8hASwhGSh5AQjpA8hIRwhOQhJIQjJA8hIRwheQgJ4QjJQ0gIR0geQkI4QvIQEsIRkoeQ\nEI6QPISEcITkISSEIyQPISEcIXkICeEIyUNICEdIHkJCOELyEBLCEZKHkBBu/qE+78vuh9KV\n1Vm1xDoICeHmHuo6H/2Ax0KyxFoICeHmHurKZcdL99H1lLlKscRaCAnh5h7qzF0eH19cplhi\nLYSEcHMP9eQHdr//6d1bO5sICeF4RPIQEsIteI50unYf8Rxp6abwBWYf6mL0ql1eS5ZYCSEh\n3IL3karufaSs3PM+0rJN4QtwZYOHkBCOkDyEhHBcIuQhJITjEiEPISEclwh5CAnheEPWQ0gI\nxyVCHkJCOB6RPISEcFwi5CEkhOMSIQ8hIRyXCHkICeG4ssFDSAgnOtRuTLOEDCEhHJcIeQgJ\n4bhEyENICMclQh5CQjjekPUQEsJxiZCHkBCORyQPISEclwh5CAnhuETIQ0gIxyVCHkJCOC4R\n8hASwhGSh5AQbvahrnfOFadhI7z8vWRT+AKzLxHK+gvt+o0Q0pJN4QvMf/n70NR0yLrL7Ahp\n0abwBea/Idv9cs3yKyEt3BS+wNJLhOqiIKSFm8IXmHuoc3d/EzYvCGnZpvAF5h7qg9sNH11d\nQUiLNoUvMPtQV496Tv/51+RbO5sICeHmH+pLef/ouiOkJZvCF+DKBg8hIRwheQgJ4QjJQ0gI\nR0geQkI4QvIQEsIRkoeQEI6QPISEcITkISSEIyQPISEcIXkICeEIyUNICEdIHkJCOELyEBLC\nEZKHkBCOkDyEhHCE5CEkhCMkDyEhHCF5CAnhCMlDSAhHSB5CQrikQnKf+WxbZkP4BmmF9NdT\n8dx0iBGSciqemw4xQlJOxXPTIUZIyql4bjrECEk5Fc9NhxghKafiuekQIyTlVDw3HWKEpJyK\n56ZDjJCUU/HcdIgRknIqnpsOMUJSTsVz0yFGSMqpeG46xAhJORXPTYcYISmn4rnpECMk5VQ8\nNx1ihKSciuemQ4yQlFPx3HSIEZJyKp6bDjFCUk7Fc9MhRkjKqXhuOsQISTkVz02HGCEpp+K5\n6RAjJOVUPDcdYoSknIrnpkOMkJRT8dx0iBGSciqemw4xQlJOxXPTIUZIyql4bjrECEk5Fc9N\nhxghKafiuekQIyTlVDw3HWKEpJyK56ZDjJCUU/HcdIgRknIqnpsOMUJSTsVz0yFGSMqpeG46\nxAhJORXPTYcYISmn4rnpECMk5VQ8Nx1ihKSciuemQ4yQlFPx3HSIEZJyKp6bDjFCUk7Fc9Mh\nRkjKqXhuOsQISTkVz02H2PxDfd6XrlVWZ9US1ggJKnMPdZ27p0KyhD1CgsrcQ1257HjpPrqe\nMlcplrBHSFCZe6gzd3l8fHGZYgl7hASVuYfaud9+Y7aEPUKCCo9Iyql4bjrEFjxHOl27j3iO\ntHRT+AKzD3UxetUuryVLmCMkqCx4H6nq3kfKyj3vIy3bFL4AVzYop+K56RAjJOVUPDcdYlwi\npJyK56ZDbHyo8/3148/jEiG7TeELTN9XdR+3xCVCdpvCFxgf6vq4+7gl3pC12xS+wOuhPu/z\nj1riEiG7TeEL/HCoL1nzuHT4z+fxiGS3KXwB/1Cfig9eQOASIULC2MuhrvfNw1F+qpuayvef\nyCVCZpvCF5gc6nP7YkPVf832/nnPjUuE7DaFLzB5H6l5bDncH1zeP++Zu8S6CAkqkxffypPZ\nZsesNroYIUFl8j5S+Kcfmkex/+UXz9lESFCZHOq6ar+ey6oPiuofZ4ZXHN6+aBfR2URIUBkf\n6mvW5eFc9v9rG7rJyrXNXav37zrFczYRElTGh7pwu/axqK7+99L3bQgpc91jV+3yT5dYFyFB\n5ccrfT54fWB47PpkPp6ziZCgMj7UwwNM8wjzYUi7e0hcIrRkU/gC40NduaJ9b/Vc/OfVg+7z\nXLk/nNzx1n0pyCVCSzaFLzA51PfLfv5znV33ec/3iJzLuERoyabwBaaH+the9VP878rvzuVy\nOJRl95LDf14uj+dsIiSo/MGhjudsIiSoEJJyKp6bDjFCUk7Fc9MhNjnU+1xxnWk8ZxMhQWV8\nqPeaC7bjOZsICSrTN2Q/er2u/7ypT5dYFyFB5cdLhD5wICSzTeELjA916QL+RdIl++BtW2+J\ndRESVKb/jKL4z7dfGLt8cCGRt8S6CAkq0y/tgl5sOIy+td2nS6yLkKCyIKQ5S6yLkKDCG7LK\nqXhuOsQISTkVz02H2PRQn8r2q7ry8x+TFL7EmggJKv6/R2r/watpSfGcTYQElfGhPrii+1fm\nB7dTLbEuQoLK9BKhevpNTeyXWBchQeX1EiFCIiTMMD7U+fCIdHn/feqWLLEuQoLKD8+RTiFX\ngQcusS5CgsrkUJeffxehuUusipCg4r+P5Mqjcok1/X1In/loQUSNKxuUU5YLImqEpJwipGQQ\nknKKkJLBP6NQThFSMghJOUVIyfjhIJ6L//+csYVLrISQoPLTQay5aNVoipCS8eNB5Es7oylC\nSsZPB/Hw/ifwWSyxDkKCys8vNuxVS6yLkKDyU0i56TWrEZ0nhAQV3pBVThFSMghJOUVIyfjl\nDVnLN2XjOU8ICSqEpJwipGRMDuI+OzX/f/7450zMWGJVhASV8UHcD98V/+JMrxGK5zwhJKhM\nv7R7/cB8iXURElTGBzF7PCLxXYRspggpGeODWLnuORLfRchsipCSMTmIxfB63Yc/im/OEqsi\nJKhMD+Kx+y5CJ+USayIkqHBlg3KKkJJBSMopQkrG9CDyg8ZspwgpGf6LDTd+0JjZFCElY3wQ\n+UFj1lOElIzpG7L8oDHbKUJKxuslQoRESJhhfBD5QWPWU4SUjB+eI3GJkNkUISVjchD5QWPG\nU4SUDP99JH7QmN0UISWDKxuUU4SUjPFBLG2v+v5piXURElReX/4WL7EuQoLK68vf4iXWRUhQ\nGR/EuizO4iXWRUhQmX5px0/ss50ipGQQknKKkJLBy9/KKUJKBiEppwgpGfeDKHrpe7zE+ggJ\nKtOQJDnFc54QElQISTlFSMkgJOUUISWDkJRThJQMQlJOEVIyCEk5RUjJeIYk+bGX4yXWR0hQ\nISTlFCElgysblFOElAxCUk4RUjIISTlFSMkgJOUUISVj/kE87/tvg1dW//lntfGcJ4QElbkH\nsc5Hr/G9/4aS8ZwnhASVuQexctnx0n10PWXvf3pzPOcJIUFl7kHM3OXx8cVliiXsERJU5h7E\nyXu279/Ajec8ISSo8IiknCKkZCx4jnTqf9Isz5H+ZkFEbfZBLEav2uVvv0NrPOcJIUFlwftI\nVfc+UlbueR/pLxZE1LiyQTlFSMkgJOUUISWDS4SUU4SUDC4RUk4RUjK4REg5RUjJ4A1Z5RQh\nJYNLhJRThJQMHpGUU4SUDC4RUk4RUjK4REg5RUjJ4BIh5RQhJYMrG5RThJQM0UGUfdvWRQgJ\nKlwipJwipGRwiZByipCSwSVCyilCSgZvyCqnCCkZXCKknCKkZPCIpJwipGRwiZByipCSwSVC\nyilCSgaXCCmnCCkZXCKknCKkZBCScoqQkkFIyilCSgYhKacIKRmEpJwipGTMv7Lh438pEc95\nQkhQmXsQD4T0xwsiarMP4iV7/48nDJYwR0hQmX8QL+8vDLJYwhohQWXBQTyMrlsVLWGMkKDC\nq3bKKUJKBiEppwgpGYSknCKkZBCScoqQkkFIyilCSgYhKacIKRmEpJwipGQQknKKkJJBSMop\nQkoGISmnCCkZhKScIqRkEJJyipCSQUjKKUJKBiEppwgpGYSknCKkZBCScoqQkkFIyilCSgYh\nKacIKRmEpJwipGQQknKKkJJBSMopQkoGISmnCCkZhKScIqRkEJJyipCSQUjKKUJKBiEppwgp\nGYSknCKkZBCScoqQkkFIyilCSgYhKacIKRmEpJwipGQQknKKkJJBSMopQkoGISmnCCkZhKSc\n+nBTn/loW1gJISmn/nxBrIWQlFOElAxCUk4RUjIISTlFSMkgJOUUISWDkJRThJQMQlJOEVIy\nCEk5RUjJICTlFCElg5CUU4SUDEJSThFSMghJOUVIySAk5RQhJYOQlFOElAxCUk4RUjIISTlF\nSMkgJOUUISWDkJRThJQMQlJOEVIyCEk5RUjJICTlFCElg5CUU4SUDEJSThFSMghJOUVIySAk\n5RQhJYOQlFOElAxCUk4RUjIISTlFSMkgJOUUISWDkJRThJQMQlJOEVIyCEk5RUjJICTlFCEl\ng5CUU4SUDEJSThFSMghJOUVIySAk5RQhJYOQlFOElAxCUk4RUjLmH57zvnStsjqrlrBGSFCZ\ne3jq3D0VkiXsERJU5h6eymXHS/fR9ZS5SrGEPUKCytzDk7nL4+OLyxRL2CMkqMw9PM799huz\nJewRElR4RFJOEVIyFjxHOl27j3iOFNGCWMvsw1OMXrXLa8kS5ggJKgveR6q695Gycs/7SNEs\niLVwZYNyipCSQUjKKUJKBpcIKacIKRlcIqScIqRkcImQcoqQksEbssopQkqG6BIhNzZzCXuE\nBBUekZRThJQMLhFSThFSMrhESDlFSMngEiHlFCElgysblFOElIxvCcl95LNN2U0RUjIWHp5D\n7lx5ki7xmTjPa0JKxrL3kYZXHN6+aEdIfzmFlSwKqXJVfbtdK3dQLBG2P2ZDhIQZFoWUue51\n79rliiXC9sdsiJAww6KQ7s/eI/guQnGe14SUjEUh7e4hrX+JUJznNSElY35I5f5wcsfmw7qK\n4BKhOM9rQkrG/JAe78w4l61/iVCc5zUhJWP24blcDoey7F5yqN52REh/OYWVfM2VDWZDhIQZ\nCEk5RUjJICTlFCElg5CUU4SUDEJSThFSMghJOUVIySAk5RQhJYOQlFOElAxCUk4RUjIISTlF\nSMkgJOWU6YJ235UC9ghJORXpgrBHSMqpSBeEPUJSTkW6IOwRknIq0gVhj5CUU5EuCHuEpJyK\ndEHYIyTlVKQLwh4hKaciXRD2CEk5FemCsEdIyqlIF4Q9QlJORbog7BGScirSBWGPkJRTkS4I\ne4SknIp0QdgjJOVUpAvCHiEppyJdEPYISTkV6YKwR0jKqUgXhD1CUk5FuiDsEZJyKtIFYY+Q\nlFORLgh7hKScinRB2CMk5VSkC8IeISmnIl0Q9ghJORXpgrBHSMqpSBeEPUJSTkW6IOwRknIq\n0gVhj5CUU5EuCHuEpJyKdEHYIyTlVKQLwh4hKaciXRD2CEk5FemCsEdIyqlIF4Q9QlJORbog\n7BGScirSBWGPkJRTkS4Ie4SknIp0QdgjJOVUpAvCHiEppyJdEPYISTkV6YKwR0jKqUgXhD1C\nUk5FuiDsEZJyKtIFYY+QlFORLgh7hKScinRB2CMk5VSkC8IeISmn/n7Bz3y0IEIQknJq0wsi\nBCEppza9IEIQknJq0wsiBCEppza9IEIQknJq0wsiBCEppza9IEIQknJq0wsiBCEppza9IEIQ\nknJq0wsiBCEppza9IEIQknJq0wsiBCEppza9IEIQknJq0wsiBCEppza9IEIQknJq0wsiBCEp\npza9IEIQknIq1gX553/mCEk59f0LYkBIyqnvXxADQlJOff+CGBCScur7F8SAkJRT378gBoSk\nnPr+BTEgJOXU9y+IASEpp75/QQwISTn1/QtiQEjKqW0vyOUPAeIPyfB6lm2f1xteMAEbCMlu\nKtLT7PsXTMD8O+K8L7uHgrI6q5YI+PQtn2bfv2AC5t4RdT76sqqQLBH06Vs+zTa9IM+kenNv\nYuWy46X76HrKXKVYIujTIz3NWLCf+v7c5u585i6Pjy8uUywR9OmbPs1YsJ/acm1z92tyg/xb\n99FN//C/U4DGzHP/55N55ucFPCIB32/Bc6TTtfvov8+RgO83++GtGD1E5rXlLgHbs+B9pKp7\nHykr9/95Hwn4frG+CAJsCiEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFC\nAgwQEmCAkAADhAQYICTAACEBBjYZ0krfvQmbJj4ntZvX0O40W19n6+z8n9vyPc7WV9r8pnde\nZMv3OFtfafOb3nmRLd/jbH2lzW9650W2fI+z9ZU2v+mdF9nyPc7WV9r8pndeZMv3OFtfafOb\n3nmRLd/jbH2lzW9650W2fI+z9ZU2v+mdF9nyPc7WV9r8pndeZMv3OFtfafOb3nmRLd/jbH2l\nzW9654FEEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBA\nSICBTYV0yF1W1e1Hgu+LXu+c2136j6vsvpBg66Lv6X4etmi/76OtK/Z9vE37nR9tXfrN9LcU\nUtXdD1lzP18E90nWbbE714vuw1yzdcW+N+qs36Jg359bV+z7eJv2Oz/auuiOH2wopIvbNQ0d\n3K79sLTeetVut+q2e3bZ5XbJ3FmydcG+t8r+DBHs+2jrin0fbVOw86Oti+74wYZCKu9fXbQ1\n7a23nrl62Hhzwp+a/z9arjHaumDfb+3e9qe6YN9HW1fs+2ibgp0fbV1zx99tKKRBH9JBtPHs\n1hZ7vUn+A9ZtXbLvV1f0p7pk3x9bV+z7aJuCnR9tXXbSdDYXUu2K9h4/7ZpnpeYbr7r72j0f\n+gRbl+x74a73J9S30S/WW1fs+2ibgp0fbV110vQ2F9Khffgv+6eNhe2mmy9gVMdztHXFvu/d\n8aYL6bl1xb6PtikJ6bF1zUlzt7WQrln7uO+aQ3urK+PH6kOZdV9Ga0J6bt1837svhmQhTbZu\nf7+PtinY+cnWFSfNYyHJVmXqbPQflNr8Vd7bTnQ8R1vvme573r4lIAtptPWe4H7vtym740d7\nrNj51sZCKib3guIez9qX2JRbHxhufde92NVv0H7fx1sfKN6Labcpu+PHmxS9kbSpkK55cR3/\nXnSP9y8eXQVvO2iOp3tQ7Pt4648/Mtv6ZJtbu+MnK0i2qnF6PFHs35Yxvcfvm8zbJ9enbjXD\nF3hGW7ff9/Gpbr/v460L7vfxNoV3fKnZ+acNhXR9vuBStfd13b9/Z6S79qAu22cxqisb+q0L\n9r0jvbLhdn+7137fR9uU3PGPravu+N6GQto9/8tY95eumb4nkD1fHM3tXyd9bl2x763hSxbB\nvj+3rtj38Tbtd360ddUd39tQSOOv1esqc7nx65jPTdbdRcjKrVvv++0RkmLfp1u33vdaesfX\n8ju+t6GQgHgREmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQUtze/wxu0U/oRjiORNwIaSM4EnEjpI3gSMSNkDaCIxGtKnPVI5Xa5d2vuatv\np9INP/u7/dt+ov//Q+6y/sd2nwrnitPf73SyCClWTQnOlY/HnMJdm/+/uuK2d522pJeQyu4v\niuajQz9zWGnfE0RIkTq67HK7ZI+Qjm7f/P/enZpmju1v27+YhnRyRX2ri2bklrlLO5OvtvvJ\nIaRIle58a9t4HKCuivz5Wz+ksvmyr/0isGx/y5d1f4uQIjUU9Axp13xtd+2+oLtdT/vih5Dc\nXfP8qvmq8HJZY79TRUiR8kI6N1/bVd3DVHHP5feQbvvmi0KXXdfY8zQRUqS8kG5Z3v6vfWjK\nD6frjyGNN3Cqcp4j/R1CilTZPcs5j+Ko3KF7waH7o5eQzv1zpJcnRrzN9He4qyN1ennVrmun\nezXBNV/fXZ7PkXJ3aF+qc8MLfbdD+2JD3r+yxyPSnyGkWHVvCu0evdzaOtq3iNoXEjrn/i8O\nz/eb+idP7TOj42MEf4OQorV/XNkwhHQcvnRr6irOp/5V7n5u97yywe26Vxi6Kxvo6O8QEmCA\nkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEsfl\nY0oAAAA2SURBVGCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgIF/a0RVWCmFVHwAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"Histogram of d.values\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the histogram of the d values\n",
    "hist(d.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case of Binomial distribution and Chi-square distribution the d-values observed are high. This is because the MLE estimates observed is seen to be really high. For other distributions we observe low values of the statistic from the ks.test function."
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}