Successive_Polynomials_Exploration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "Fit a series of successive polynomial models to the data: i.e. first fit\n",
    "a constant, then fit a linear model, then fit a quadratic model, then a cubic model, etc. What\n",
    "is the degree of the first degree of the first model that you would accept, at the 1% confidence\n",
    "level? What is the degree of the first model that you would accept at the 1% confidence level? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xi:  [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20]\n",
      "Y:  [141 144 148 150 158 161 166 170 175 181 189 194 196 206 218 229 234 242\n",
      " 247 257]\n"
     ]
    }
   ],
   "source": [
    "xi = np.array([i for i in range(1, 21)])\n",
    "print(\"xi: \", xi)\n",
    "\n",
    "Y = np.array([141, 144, 148, 150, 158, 161, 166, 170, 175, 181, 189, 194, 196, 206, 218, 229, 234, 242, 247, 257])\n",
    "print(\"Y: \", Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def orthog_poly(n, xi):\n",
    "    # Define the first two orthogonal polynomials\n",
    "    psi = []\n",
    "    phi = []\n",
    "    psi.append(lambda t: np.ones_like(t))\n",
    "    psi.append(lambda t: t - np.mean(t))\n",
    "\n",
    "    for i in range(2, n+1):\n",
    "        a = np.sum(xi*(psi[i-1](xi))**2)/np.sum((psi[i-1](xi))**2)\n",
    "        b = np.sum(psi[i-1](xi)**2)/np.sum(psi[i-2](xi)**2)\n",
    "        \n",
    "        psi.append(lambda t, i=i, a=a, b=b: (t - a)*psi[i-1](t) - b*psi[i-2](t))\n",
    "    \n",
    "    for i in range(n + 1):\n",
    "        phi.append(lambda t, i=i: psi[i](t)/np.sqrt(np.sum(psi[i](xi)**2)))    \n",
    "    \n",
    "    return phi\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poly_size:  1\n",
      "p_value:  1.1102230246251565e-16\n",
      "Reject null hypothesis\n",
      "\n",
      "poly_size:  2\n",
      "p_value:  1.1007069034008055e-08\n",
      "Reject null hypothesis\n",
      "\n",
      "poly_size:  3\n",
      "p_value:  0.6458808944859353\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  4\n",
      "p_value:  0.07515056187201463\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  5\n",
      "p_value:  0.1597879250036146\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  6\n",
      "p_value:  0.36960429118117355\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  7\n",
      "p_value:  0.08863805302818262\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  8\n",
      "p_value:  0.12385260108166918\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  9\n",
      "p_value:  0.6383626921996918\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  10\n",
      "p_value:  0.07778119221002755\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  11\n",
      "p_value:  0.31868242856920914\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  12\n",
      "p_value:  0.05072620721736554\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  13\n",
      "p_value:  0.12413741051150229\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  14\n",
      "p_value:  0.1225909170386924\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  15\n",
      "p_value:  0.04467298247549267\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  16\n",
      "p_value:  0.38599834816926537\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  17\n",
      "p_value:  0.5461356465261572\n",
      "Fail to reject null hypothesis\n",
      "\n",
      "poly_size:  18\n",
      "p_value:  0.0022459734453023072\n",
      "Reject null hypothesis\n",
      "\n",
      "poly_size:  19\n",
      "p_value:  1.1618387696366028e-08\n",
      "Reject null hypothesis\n",
      "\n"
     ]
    }
   ],
   "source": [
    "p_value = .01\n",
    "\n",
    "# Fit a series of successive polynomial models to the data\n",
    "for poly_size in range(1, len(Y)):\n",
    "    orthog_poly_basis = orthog_poly(poly_size, xi)\n",
    "    X_orth =np.array([q(xi) for q in orthog_poly(poly_size, xi)]).T\n",
    "    beta_hat_orth = X_orth.T @ Y\n",
    "    \n",
    "    RSS = np.sum((Y - X_orth @ beta_hat_orth)**2)\n",
    "    #print(\"RSS: \", RSS)\n",
    "    #print(\"beta_hat_orth: \", beta_hat_orth)\n",
    "    \n",
    "    # under null hypothesis that the leading coefficient is zero\n",
    "    A = np.zeros((1, poly_size + 1))\n",
    "    A[:, -1] = 1\n",
    "\n",
    "    c = 0\n",
    "    beta_h0 = beta_hat_orth + A.T@ np.linalg.inv(A @ A.T) @ (c -  A @ beta_hat_orth)\n",
    "    # print(\"beta_h0: \", beta_h0)\n",
    "    # RSS_h0 = np.sum((Y -  X_orth @ beta_h0 )**2)\n",
    "    RSS_h0 = (Y -  X_orth @ beta_h0 ).T @ (Y -  X_orth @ beta_h0 )\n",
    "    #print(\"RSS_h0: \", RSS_h0)\n",
    "    \n",
    "    # calculate f stat\n",
    "    f_stat_crit = np.abs(RSS_h0 - RSS)/RSS * (len(Y) - poly_size)\n",
    "\n",
    "    #print(\"diff: \", RSS_h0 - RSS)\n",
    "    #print(\"f_stat: \", f_stat_crit)\n",
    "    print(\"poly_size: \", poly_size)\n",
    "    \n",
    "    p_compare =  1 - scipy.stats.f.cdf(f_stat_crit, 1, len(Y) - poly_size)\n",
    "    print(\"p_value: \", p_compare)\n",
    "    if p_value > p_compare:\n",
    "        print(\"Reject null hypothesis\\n\")\n",
    "    else:\n",
    "        print(\"Fail to reject null hypothesis\\n\")\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The degree of the first model I would accept at the 1% confidence level is 2. See above. \n",
    "\n",
    "This is an odd question because it asks whether the improvement over the previous model is statistically significant. So, obviously a large improvement would be made from the linear model to the quadratic. It is odd that this also happens at degree = 19. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}