2-Probability and Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "40b8b45a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import exists\n",
    "import requests\n",
    "from typing import *\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import scipy.stats as stats\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm.notebook import trange, tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c3591ce",
   "metadata": {},
   "source": [
    "# Notebook 2: Probability and Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "916d1cfa",
   "metadata": {},
   "source": [
    "## Problem 1: A Strange Regression Model\n",
    "\n",
    "Suppose we are fitting a regression model to a dataset $(x_i, y_i)_{1 \\leq i \\leq N}$\n",
    "\\begin{align*}\n",
    "p(y^i|x^i; \\theta) & = \\mathcal{N}(f_1(x_1^i) \\theta_1 + f_2(x_2^i) \\theta_2 + \\theta_3, 1) \\\\\n",
    "p(y |x; \\theta) & = \\prod_{i=1}^N p(y^i|x^i; \\theta)\n",
    "\\end{align*}\n",
    "where\n",
    "1. the inputs $x^i \\in \\mathbb{R}^2$ are 2-D vectors\n",
    "2. $f_1, f_2, f_3: \\mathbb{R} \\rightarrow \\mathbb{R}$ are arbitrary functions.\n",
    "3. $\\theta = (\\theta_1, \\theta_2, \\theta_3)$ is a vector of weights"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e48dd24a",
   "metadata": {},
   "source": [
    "### Problem 1a\n",
    "\n",
    "Implement the conditional density given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "62d050f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scout resources to verify correctness of the implementation.\n",
    "\n",
    "x1 = [2., 2.]\n",
    "x2 = [3., 4.]\n",
    "x3 = [5., 6.]\n",
    "\n",
    "x = np.array([x1, x2, x3])\n",
    "Y = np.array([7., 13., 30.])\n",
    "try_theta = np.array([.3, .67, .25])\n",
    "actual_thetas = [1.0, 0.5, 2.]\n",
    "\n",
    "def f1(x: float) -> float:\n",
    "    return x**2\n",
    "\n",
    "def f2(x: float) -> float:\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e74537c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def model1_density(f1: Callable[[float], float], f2: Callable[[float], float], theta: np.ndarray, x: np.ndarray) -> Callable[[np.ndarray], float]:\n",
    "    # f1 is function f_1\n",
    "    # f2 is function f_2\n",
    "    # theta is some choice of weights\n",
    "    # x is an input\n",
    "    # y: length N vector\n",
    "\n",
    "    x_new = np.concatenate([x, np.ones((x.shape[0], 1))], axis=1)\n",
    "    x_new[:, 0] = f1(x_new[:, 0])\n",
    "    x_new[:, 1] = f2(x_new[:, 1])\n",
    "\n",
    "    normals = [ stats.norm(loc= theta @ x_new[i], scale=1) for i in range(x_new.shape[0]) ]\n",
    "\n",
    "    #conditional_density = normals[0].pdf(Y[0]) * normals[1].pdf(Y[1]) * normals[2].pdf(Y[2])\n",
    "\n",
    "    return lambda y: np.prod([normals[i].pdf(y[i]) for i in range(len(y))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "720d2377",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.864009903343106e-90\n"
     ]
    }
   ],
   "source": [
    "density = model1_density(f1, f2, try_theta, x)\n",
    "\n",
    "print(density(Y))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "068d0eee",
   "metadata": {},
   "source": [
    "### Problem 1b\n",
    "\n",
    "Write a function that solves for the weights by finding the approximate minimum of the conditional density, i.e., solve\n",
    "$$\n",
    "\\operatorname{argmin}_{\\theta} -p(y | x; \\theta)\n",
    "$$\n",
    "1. The input $X^{N \\times 2}$ is a $N \\times 2$ matrix where row $i$ of the matrix contains the input $x^i$.\n",
    "2. The input $Y^N$ is a length $N$ vector where each entry contains the corresponding value $y^i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "289ecd8b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f1618262be940cab952adcbcd0d8234",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c038564280a48af9fde88738b13a57f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Thetas obatained from grid search: [0.84848485 1.45454545 0.        ]\n",
      "Probability density based on these thetas: 0.04483163226933461\n"
     ]
    }
   ],
   "source": [
    "def solve_for_weights(f1: Callable[[float], float], f2: Callable[[float], float], X: np.ndarray, Y: np.ndarray) -> np.ndarray:\n",
    "    \n",
    "    thetaspace = np.linspace(0., 12.0, 100)\n",
    "\n",
    "    best_density = 0.\n",
    "    best_theta = np.array([1.5, .25, 0])\n",
    "    \n",
    "    for theta3 in tqdm(thetaspace):\n",
    "        theta = np.array([best_theta[0], best_theta[1], theta3])\n",
    "        density_f = model1_density(f1, f2, theta, x)\n",
    "        density = density_f(Y)\n",
    "        if density > best_density:\n",
    "            best_density = density\n",
    "            best_theta[2] = theta3\n",
    "\n",
    "    for theta2 in tqdm(thetaspace):\n",
    "        for theta1 in thetaspace:\n",
    "            theta = np.array([theta1, theta2, best_theta[2]])\n",
    "            density_f = model1_density(f1, f2, theta, x)\n",
    "            density = density_f(Y)\n",
    "            if density > best_density:\n",
    "                best_density = density\n",
    "                best_theta[0] = theta1\n",
    "                best_theta[1] = theta2\n",
    "\n",
    "    return best_theta\n",
    "\n",
    "grid_thetas = solve_for_weights(f1, f2, x, Y)\n",
    "\n",
    "print(\"Thetas obatained from grid search: \" + str(grid_thetas))\n",
    "\n",
    "density_f = model1_density(f1, f2, grid_thetas, x)\n",
    "density = density_f(Y)\n",
    "print(\"Probability density based on these thetas: \" + str(density))\n",
    "\n",
    "#can try plotting the density or set functions as identity functions: x comonents themselves"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb5a24b",
   "metadata": {},
   "source": [
    "### Problem 1c\n",
    "\n",
    "Transform the probabilistic model into a system that can be solved with linear algebra.\n",
    "```ts\n",
    "def linear_algebra_solve(X, Y):\n",
    "    X_p = np.concatenate([X, np.ones((Y.size[0], 1))], axis=1)\n",
    "    return np.linalg.solve(X_p.transpose() @ X_p, X_p.transpose() @ Y)\n",
    "```\n",
    "That is, come up with new matrix $X'$ and vector $Y'$ such that it can be input into `linear_algebra_solve`ts to produce the same set of weights.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "38d3e4c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Thetas obtained from linear algebra: [1.  0.5 2. ]\n",
      "Probability density based on these thetas: 0.06349363593424098\n"
     ]
    }
   ],
   "source": [
    "def transform(f1: Callable[[float], float], f2: Callable[[float], float], X: np.ndarray, Y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:\n",
    "    \n",
    "    # X is a matrix of inputs\n",
    "    # Y is a vector of outputs\n",
    "\n",
    "    # return a new X and Y\n",
    "\n",
    "    new_X = np.array([f1(x[:, 0]), f2(x[:, 1])]).T\n",
    "    #new_Y = Y\n",
    "\n",
    "    return new_X, Y\n",
    "\n",
    "def linear_algebra_solve(X, Y):\n",
    "    X_p = np.concatenate([X, np.ones((X.shape[0], 1))], axis=1)\n",
    "    return np.linalg.solve(X_p.transpose() @ X_p, X_p.transpose() @ Y)\n",
    "\n",
    "las_X, las_Y = transform(f1, f2, x, Y)\n",
    "linalg_thetas = linear_algebra_solve(las_X, las_Y)\n",
    "\n",
    "print(\"Thetas obtained from linear algebra: \" + str(linalg_thetas))\n",
    "\n",
    "linalg_density_f = model1_density(f1, f2, linalg_thetas, x)\n",
    "\n",
    "print(\"Probability density based on these thetas: \" + str(linalg_density_f(Y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4526ee4",
   "metadata": {},
   "source": [
    "## Problem 2: Mixture Model\n",
    "\n",
    "Suppose that we are performing regression in a setting with the following density\n",
    "$$\n",
    "p(y|x; \\theta, \\beta) = \\frac{1}{2}p_1(y|x; \\theta) + \\frac{1}{2}p_2(y|x; \\beta)\n",
    "$$\n",
    "where both $p_1(y|x; \\theta)$ and $p_2(y|x; \\theta)$ are linear regressions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55784262",
   "metadata": {},
   "source": [
    "### Problem 2a\n",
    "\n",
    "Write a function that generates synthetic data according to the following probabilistic model.\n",
    "\\begin{align*}\n",
    "z^i & \\sim Bern(0.5) \\\\\n",
    "y^i & \\sim \\begin{cases}\n",
    "\\mathcal{N}(\\theta x^i, 1) & \\mbox{when $z^i = 0$} \\\\\n",
    "\\mathcal{N}(\\beta x^i, 1) & \\mbox{when $z^i = 1$}\n",
    "\\end{cases}\n",
    "\\end{align*}\n",
    "In this example, we are assuming that the inputs are a single real number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0c337a40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcJElEQVR4nO3dfZAc9X3n8feHZQ8WsL1SsWBpJZBypwjzEKPUgknkSnwCTnIgSKaOGHwmMqZO5yr7/FBEvlXsuuCcE5RSznHOia9OsTFy4mAIpoQCvlNkKTaxzINXrM4ghA5snrQSaHlYHuw11sP3/uhePNqZnt0ZzUzPw+dVtbU73T3d3xFUf+f39G1FBGZmZoWOyzsAMzNrPk4OZmZWxMnBzMyKODmYmVkRJwczMyvi5GBmZkWcHKxjSXpK0iU1OtcZkl6X1FWL85nlzcnBmoqkd0v6gaRXJL0kabukC2pw3lskfb4WMabnOyqxRMQzEXFKRByu8DwfknQ4TSyvS3pS0tck/WoF56jpZ8v7OtYcnBysaUh6K3A38CVgJtAPfA54I8+4GuC+iDgFeBtwCTAO7JB0br5hWUeLCP/4pyl+gAFgLGPfCcBLwHkF204juZH2Ae8B9gI3AAeA/cB16XGrgIPAL4DXgX9Mtz8F/AHwI+AV4DbgxILzXw7sBMaAHwC/lm7/W+BIeu3XgU8D84AAjk+PmQl8DdgHvAxszPhcHwK+X2L73cAdBa//AXgujfNe4JwpPtsg8GPgNeBR4H0F5/o3wPfSc70A3Faw7yxgS/pvvQf4vXLX8U/7/uQegH/8M/EDvBV4EdgAvBeYMWn/l4E/K3j9iYKb4XuAQ8AfA93A7wA/mzgHcAvw+Unnewp4EJid3sx3Ax9J9/16mmTeBXQBK9PjTyh47yUF55qcHO5Jk82MNJ7fzvjMWcnhw8Dzk16/hSRJfhHYWbCv1Ge7Kv1cxwHvB34KzEr33Qp8Jt13IvDudPvJwLPAdcDx6b/BCwWJqOg6/mnfH3crWdOIiFeBd5PcZP8GGJW0SdLp6SEbgA9Imvj/9lqSb/ETDgJ/HBEHI+LbJN9wF05x2f8REfsi4iXgH4Hz0+3/EfhfEfFARByOiA0k3VsXTfU5JM0iSW4fiYiX03i+N9X7JtlHkrAAiIibI+K1iHgDuBF4p6S3Zb05Iv4h/VxHIuI24HHgwnT3QeBMYHZE/Dwivp9uvxx4KiK+FhGHIuIh4FvAv68wdmsDTg7WVCJid0R8KCLmAOeSfPv9YrrvAZJvwL8t6SyS7pFNBW9/MSIOFbz+GXDKFJd8LuP4M4EbJI1N/ABz03imMhd4KSJensaxWfpJunaQ1CVpraQfS3qVpNUCcGrWmyX9vqSdBbGfW3D8pwEBD0raJenD6fYzgXdN+sz/AXj7MXwOa1HH5x2AWZaIeEzSLcB/Kti8AfggyU39joj4+XRPV+HlnwX+JCL+pIrzPQvMlNQbEWMVXnfC+4B/Sf/+ALCcZLD6KZKB65dJbvBFsUg6k6TldTHJYPdhSTsnjo+I50haRkh6N/AdSfemcX8vIi7NiMklnDuIWw7WNCSdJekGSXPS13OBa4D7Cw77W5Ib5weBr1dw+ueBX6ng+L8BPiLpXUqcLOkySW+Z6nwRsR/438CXJc2Q1C3pt6a6YNpCmC/pSyRjKJ9Ld72FpEvrReAk4E+n+Gwnk9zIR9PzXkfScpi4zlUT/8YkSSaAwySD4L8q6do05m5JF0h6x1Sf2dqPk4M1k9dIBoAfkPRTkqTwCMkMJAAiYi/wEMkN7V9KnSTDV4Gz0+6SjVMdHBFDJN+u/4rkBvoEyeDxhJuAz6bn+4MSp7iWpG//MZKB7U+WudxvSHodeBX4LsnA/AUR8XC6/+vA08AIycyj+ye9/6jPFhGPAv8duI/khn4esL3g+AtI/o1fJ+mW+0REPBkRrwH/DriaZMzjOeDPSAbBi65T5vNYG1CEW4rWWiTdDOyLiM/mHYtZu/KYg7UUSfOAK4FFOYdi1tbcrWQtQ9J/I+lmWhcRT+Ydj1k7c7eSmZkVccvBzMyKtMWYw6mnnhrz5s3LOwwzs5ayY8eOFyKir9S+tkgO8+bNY2hoKO8wzMxaiqSns/bl3q2ULvwZlnR3+nqmpC2SHk9/z8g7RjOzTpN7ciCprLm74PUgsDUiFgBb09dmZtZAuSaHdAn/ZcBXCjYvJ6mfQ/p7RYPDMjPreHm3HL5IUiHySMG209PaNBM1ak4r9UZJqyQNSRoaHR2te6BmZp0kt+Qg6XLgQETsqOb9EbE+IgYiYqCvr+Rgu5mZVSnP2UqLgSsk/Q7J06jeKunvgOclzYqI/elDUw7kGKOZWd1tHB5h3eY97BsbZ3ZvD6uXLmTFov5cY8qt5RARayJiTkTMI6kCuS0iPkhSJXJlethK4K6cQjQzq7uNwyOsufNhRsbGCWBkbJw1dz7MxuGRXOPKe8yhlLXApZIeBy5NX5uZtaV1m/cwfvDwUdvGDx5m3eY9OUWUaIpFcBHxXZI69kTEiyRPsDIza3v7xsYr2t4ozdhyMDPrGLN7eyra3ihODmZmOVq9dCE93V1Hbevp7mL10oU5RZRoim4lM7NONTErqdlmKzk5mJnlbMWi/tyTwWTuVjIzsyJODmZmVsTdSmZmVWrGlc214uRgZlaFiZXNEwvYJlY2A22RINytZGZWhWZd2VwrTg5mZlVo1pXNteLkYGZWhWZd2VwrHnMwM5tCqYHn1UsXHjXmAM2xsrlW3HIwMysjq6Q2wE1Xnkd/bw8C+nt7uOnK89piMBrccjAzK6vcwPP2wSVtkwwmc8vBzKyMdh94zuKWg5m1tHovRJvd28NIiUTQLgPPWdxyMLOW1YhHbDZrSe16c3Iws5bViIVoKxb1t/XAcxZ3K5lZy2rUeEAzltSuN7cczKxltftCtDw5OZhZy+rU8YBGcLeSmbWsZn3EZjtwcjCzltaJ4wGN4G4lMzMrkltykHSipAcl/V9JuyR9Lt0+U9IWSY+nv2fkFaOZWafKs+XwBrAkIt4JnA8sk3QRMAhsjYgFwNb0tZmZNVBuySESr6cvu9OfAJYDG9LtG4AVjY/OzKyz5TrmIKlL0k7gALAlIh4ATo+I/QDp79NyDNHMrCPlmhwi4nBEnA/MAS6UdO503ytplaQhSUOjo6N1i9HMrBM1xWyliBgDvgssA56XNAsg/X0g4z3rI2IgIgb6+voaFaqZWUfIc7ZSn6Te9O8e4BLgMWATsDI9bCVwVy4Bmtkx2zg8wuK125g/eA+L126rabVUq688F8HNAjZI6iJJUrdHxN2S7gNul3Q98AxwVY4xmlmVJsppT1RNLXy8phetNb/ckkNE/AhYVGL7i8DFjY/IzGqpXDltJ4fm1xRjDmbWfjr18ZrtwsnBzOrC5bRbm5ODmdWFy2m3NldlNbO6cDnt1ubkYGZ143LarcvJwcyO2cbhEbcQ2oyTg5kdE69naE9ODmZtrt7f6r2eoT05OZi1sUZ8q/d6hvbkqaxmbazct/pa8XqG9uTkYNbGGvGt3usZ2pOTg1kba8S3+hWL+rnpyvPo7+1BQH9vDzddeZ7HG1qcxxzM2tjqpQuPGnOA+nyr93qG9uPkYNbGvErZquXkYNbm/K3equExBzMzK+LkYGZmRZwczMysiMcczHJUaWkLF7izRnFyMMtJpaUtXODOGsndSmY5qbS0RSNKYZhNcHIwy0mlpS1c4M4aycnBLCeVlrZwgTtrJCcHs5xUWrDOBe6skTwgbZaTSktbuBSGNZIiIp8LS3OBrwNvB44A6yPiLyXNBG4D5gFPAb8XES+XO9fAwEAMDQ3VN2AzszYjaUdEDJTal2e30iHghoh4B3AR8FFJZwODwNaIWABsTV+bmVkD5ZYcImJ/RDyU/v0asBvoB5YDG9LDNgArcgnQzKyDNcWAtKR5wCLgAeD0iNgPSQIBTst4zypJQ5KGRkdHGxarmVknyD05SDoF+BbwyYh4dbrvi4j1ETEQEQN9fX31C9DMrAPlmhwkdZMkhm9ExJ3p5uclzUr3zwIO5BWfmVmnyi05SBLwVWB3RHyhYNcmYGX690rgrkbHZmbW6fJc57AYuBZ4WNLOdNsfAmuB2yVdDzwDXJVPeGZmnSu35BAR3weUsfviRsZiZmZHy31A2szMmo+Tg5mZFXFyMDOzIk4OZmZWxMnBzMyKODmYmVkRJwczMyvih/1Yx9g4POIH5ZhNk5ODdYSNwyOsufNhxg8eBmBkbJw1dz4M4ARhVoK7lawjrNu8583EMGH84GHWbd6TU0Rmzc3JwTrCvrHxirabdTonB+sIs3t7Ktpu1umcHKypbRweYfHabcwfvIfFa7excXikqvOsXrqQnu6uo7b1dHexeunCWoRp1nY8IG1Nq5aDyBPHe7aS2fQ4OVjTKjeIXM1NfcWificDs2lyt5I1LQ8im+XHycGalgeRzfLj5GBNy4PIZvnxmIM1LQ8im+XHycGamgeRzfLhbiUzMyvi5GBmZkWcHMzMrIiTg5mZFXFyMDOzIrkmB0k3Szog6ZGCbTMlbZH0ePp7Rp4xmpl1orxbDrcAyyZtGwS2RsQCYGv62szMGijX5BAR9wIvTdq8HNiQ/r0BWNHImMzMLP+WQymnR8R+gPT3aaUOkrRK0pCkodHR0YYGaGbW7poxOUxLRKyPiIGIGOjr68s7HDOztjJlcpD0sQYPCj8vaVZ67VnAgQZe28zMmF7L4e3ADyXdLmmZJNU5pk3AyvTvlcBddb6emZlNMmVyiIjPAguArwIfAh6X9KeS/vWxXlzSrcB9wEJJeyVdD6wFLpX0OHBp+trMzBpoWlVZIyIkPQc8BxwCZgB3SNoSEZ+u9uIRcU3GrourPaeZmR27KZODpI+TdO+8AHwFWB0RByUdBzwOVJ0czMysOU2n5XAqcGVEPF24MSKOSLq8PmHZZBuHR/zQGzNrmCmTQ0T81zL7dtc2HCtl4/AIa+58mPGDhwEYGRtnzZ0PAzhBmFldtOw6h06ybvOeNxPDhPGDh1m3eU9OEZlZu3NyaAH7xsYr2m5mdqycHFrA7N6eirabmR0rJ4cWsHrpQnq6u47a1tPdxeqlC3OKyMza3bTWOVi+JgadPVvJzBrFyaFFrFjU72RgZg3jbiUzMyvilkNOsha1ebGbmTUDJ4ccZC1qG3r6Jb61Y8SL3cwsd+5WykHWorZbH3jWi93MrCk4OeQga/Ha4YiKjjczqxcnhxxkLV7ryniOkhe7mVmjOTnkIGtR2zXvmuvFbmbWFDwgXUdZM4/KLWobOHOmZyuZWe4UGf3crWRgYCCGhobyDuMok2ckQdIKuOnK83yzN7OmIGlHRAyU2udupTpxmW0za2VODnXiMttm1sqcHOrEZbbNrJU5OdSJy2ybWSvzbKUKVFL3yGW2zayVOTlMU1Y9JMiue+Qy22bWqjo6OVTSEig3+8gJwMzaTdOOOUhaJmmPpCckDdb6/BMtgZGxcYJftgQ2Do+UPN6zj8yskzRlcpDUBfw18F7gbOAaSWfX8hqVrkPw7CMz6yTN2q10IfBERPwEQNI3geXAo7W6QLmWQKnuptVLF5Zc8ezZR2bWjpqy5QD0A88WvN6bbnuTpFWShiQNjY6OVnyBrG/8b+vpLtndBHDTlefR39uDgP7eHpfCMLO21awth1K1q48qAhUR64H1kNRWqvQCWS0Biczupu2DS5wMzKwjNGvLYS8wt+D1HGBfLS+wYlF/yZbA2M8Oljy+3MDzxuERFq/dxvzBe1i8dlvmoLaZWato1pbDD4EFkuYDI8DVwAdqfZFS6xDWbd7DSIlEkNUNVc36BzOzZteULYeIOAR8DNgM7AZuj4hdjbh2pWUvXH3VzNpRs7YciIhvA99u9HUrLXuR9/qHShbymZlNV9MmhzxVUvZidm9PRd1QteQuLTOrl6bsVmoleVZfrbZLywPoZjYVtxyOUZ7VV6vp0nJrw8ymw8mhBvKqvlquSytrLMIFBM1sOtyt1MKyurT+7Vl9mUUF8x5AN7PW4OTQwrIW8v3zY6OZrQMXEDSz6XC3Uosr1aX1qdt2ljx239g4f/H+82tWQNDTaM3al5NDGyo3FlHNAHqpJAB4YNusjSmi4pp1TWdgYCCGhobyDqNpTJ6RBEnroJoqslnnOrH7OF4uUYeqv7eH7YNLqg/ezBpG0o6IGCi1zy2HNlTL6bVZs5smb5vggW2z9uDk0KaqmV5bqvuo0pu9B7bN2oNnKxmQ/Uzt3pO6Sx7f29Od28pwM6s/txwMyO4+OuH44+jp7ioac7jxinPefF8tZit55pNZc3FyMCB7rOCV8YP8xfvPz7xx1+IG7pIeZs3HycGAqae/1mL8IuscLulh1nycHAzIfqZ2tYvjsloCUNwVVeuSHu6iMjt2Tg4GNGb6642bdvHGoSNFSaP3pO6SaybKzXzKSgDuojKrDScHe1OtqstmfeMfGy9OAOUGvbNaLeUSgLuozGrDU1mt5ipd6/DK+MGSBQSrGaNw1Vmz2nDLwWoua/wiq+RGpYPe5RJAno9tNWsnbjlYzWWVEv+j3z2nJgvnypUdr+axrX5sqlkxtxysLsq1BI510LvczKpKB9YrnVnlcQvrFK7Kai2pVtNVF6/dVrIbqren+6iZVVB9ZVuzZuWqrNZ28ppZNdWsJ6+xsHaRy5iDpKsk7ZJ0RNLApH1rJD0haY+kpXnEZ52j0oHqfWPjmWMUWcULPYZhrSivAelHgCuBews3SjobuBo4B1gGfFlSV/HbzWojawB7RkY12rf1dGcmgHJTbM1aTS7dShGxG0DS5F3LgW9GxBvAk5KeAC4E7mtshNYpsgawgZKD3hJVrbFoRHeTu7SslpptzKEfuL/g9d50WxFJq4BVAGeccUb9I7O2VcnMqk/dtrPkceXWWEy0Niop6VHpjd5lQ6zW6tatJOk7kh4p8bO83NtKbCs5nSoi1kfEQEQM9PX11SZoswIrFvWzfXAJT669jO2DS1ixqL+qNRblWhulVDN24S4tq7W6JYeIuCQizi3xc1eZt+0F5ha8ngPsq1eMZpUqt8gua/HfWIlV4ZA9U6qaG73LhlitNVu30ibg7yV9AZgNLAAezDcks1+aapFdqS6qdZv3VFTSo5obfS3LhnjswiCn5CDpfcCXgD7gHkk7I2JpROySdDvwKHAI+GhEHC53LrNGq3SNRaXPyqjmRl/N8zhKJQHAYxcGeIW0WUNU8m188uAyTG91di2ukVUcsb+3h+2DS6b7ca1FlFsh7eRg1oTq3bWTVTYki4An115Ws+u766o5uHyGWYupVXmQLJUOVNey5Lmn3bYGl+w260BZN/venu6alFUvx9NuW4OTg1kHypqSe+MV51T0VL5qeNpta3C3klkHms6U3Hrx0/pag5ODWYeq97hGlmqm3VrjOTmYWUNV+rQ+y4eTg5lNW62moObVarHpc3Iws2nxFNTO4tlKZjYtnoLaWdxyMLNpacQUVK+cbh5ODmY2LfWegjpVt1U1D0Byoqmek4OZTUu9p6BO1W1VyXhHuUQzcS0njfKcHMxsWuo9BbVct1W5xFHq+lnH37hpF28cOpLboHortWacHMxs2uo5BbVct1Wl4x1Z28fGi8uRF7ZO6nnjbrXZXp6tZGZNodwjWMs9u7uS7VkmbtSVPLe7Uq0228vJwcyaQtYzuFcs6i+bOErJOn7GSd0lj++S6n7jbrWCg+5WMrOmkdVtVel4R9bxQMlB9cmJYUItb9ytVnDQycHMWkKl4x3ljp+cNNZt3lPTG3epgedWKzjox4SaWcer9rndlZ4LmmsarR8TamZWRi2n6ZYbeN4+uKQpZyaV4uRgZkbtpuk2auC53msmPFvJzKyGajW9tpyJrqt6Tr11cjAzq6FKp91WoxFrJtytZGZWQ4140l0juq5ySQ6S1gG/C/wC+DFwXUSMpfvWANcDh4GPR8TmPGI0M6tWvZ9014g1E3l1K20Bzo2IXwP+H7AGQNLZwNXAOcAy4MuSujLPYmbW5jYOj7B47TbmD97D4rXb2Dg80pCuq1ySQ0T8U0QcSl/eD8xJ/14OfDMi3oiIJ4EngAvziNHMLG9ZA89AZqmRWmmGMYcPA7elf/eTJIsJe9NtRSStAlYBnHHGGfWMz8wsF3mumahbcpD0HeDtJXZ9JiLuSo/5DHAI+MbE20ocX3IJd0SsB9ZDskL6mAM2M2syeRbrq1tyiIhLyu2XtBK4HLg4flnDYy8wt+CwOcC++kRoZtbc8izWl8uYg6RlwH8BroiInxXs2gRcLekESfOBBcCDecRoZpa3Rgw8Z8lrzOGvgBOALZIA7o+Ij0TELkm3A4+SdDd9NCJK19I1M2tzjVgzkcVVWc3MOlS5qqwun2FmZkWcHMzMrIiTg5mZFXFyMDOzIk4OZmZWpC1mK0kaBZ4+hlOcCrxQo3BaQad9XvBn7hT+zJU5MyL6Su1oi+RwrCQNZU3naked9nnBn7lT+DPXjruVzMysiJODmZkVcXJIrM87gAbrtM8L/sydwp+5RjzmYGZmRdxyMDOzIk4OZmZWpKOTg6RlkvZIekLSYN7x1JukuZL+WdJuSbskfSLvmBpFUpekYUl35x1LI0jqlXSHpMfS/96/kXdM9STpU+n/049IulXSiXnHVA+SbpZ0QNIjBdtmStoi6fH094xaXKtjk4OkLuCvgfcCZwPXSDo736jq7hBwQ0S8A7gI+GgHfOYJnwB25x1EA/0l8H8i4izgnbTxZ5fUD3wcGIiIc4Eu4Op8o6qbW4Blk7YNAlsjYgGwNX19zDo2OQAXAk9ExE8i4hfAN4HlOcdUVxGxPyIeSv9+jeSGUf+nhuRM0hzgMuArecfSCJLeCvwW8FWAiPhFRIzlGlT9HQ/0SDoeOIk2fbxwRNwLvDRp83JgQ/r3BmBFLa7VycmhH3i24PVeOuBGOUHSPGAR8EDOoTTCF4FPA0dyjqNRfgUYBb6WdqV9RdLJeQdVLxExAvw58AywH3glIv4p36ga6vSI2A/JF0DgtFqctJOTg0ps64h5vZJOAb4FfDIiXs07nnqSdDlwICJ25B1LAx0P/DrwPyNiEfBTatTV0IzSPvblwHxgNnCypA/mG1Xr6+TksBeYW/B6Dm3aFC0kqZskMXwjIu7MO54GWAxcIekpkq7DJZL+Lt+Q6m4vsDciJlqFd5Aki3Z1CfBkRIxGxEHgTuA3c46pkZ6XNAsg/X2gFift5OTwQ2CBpPmS/hXJANamnGOqK0ki6YfeHRFfyDueRoiINRExJyLmkfw33hYRbf2tMiKeA56VtDDddDHwaI4h1dszwEWSTkr/H7+YNh6AL2ETsDL9eyVwVy1OenwtTtKKIuKQpI8Bm0lmN9wcEbtyDqveFgPXAg9L2plu+8OI+HZ+IVmd/GfgG+kXn58A1+UcT91ExAOS7gAeIpmRN0ybltGQdCvwHuBUSXuBPwLWArdLup4kUV5Vk2u5fIaZmU3Wyd1KZmaWwcnBzMyKODmYmVkRJwczMyvi5GBmZkWcHMzMrIiTg5mZFXFyMKsDSRdI+pGkEyWdnD5r4Ny84zKbLi+CM6sTSZ8HTgR6SGod3ZRzSGbT5uRgVidp6YofAj8HfjMiDuccktm0uVvJrH5mAqcAbyFpQZi1DLcczOpE0iaSMuHzgVkR8bGcQzKbto6tympWT5J+HzgUEX+fPq/8B5KWRMS2vGMzmw63HMzMrIjHHMzMrIiTg5mZFXFyMDOzIk4OZmZWxMnBzMyKODmYmVkRJwczMyvy/wFsiUSYiDyG/QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def generate_dataset(theta: float, beta: float, x: np.ndarray) -> np.ndarray:\n",
    "    # X is a length N array\n",
    "    # Return N y's as a np.\n",
    "    \n",
    "    bern = sp.stats.bernoulli(0.5)\n",
    "    z = bern.rvs(size=50)\n",
    "\n",
    "    z0 = lambda x: stats.norm(loc=theta * x, scale=1)\n",
    "    z1 = lambda x: stats.norm(loc=beta * x, scale=1)\n",
    "\n",
    "    y = np.array([z0(x[i]).rvs() if z[i] == 0 else z1(x[i]).rvs() for i in range(len(x))])\n",
    "\n",
    "    #for i in range(len(z)):\n",
    "    #    if z[i] == 0:\n",
    "    #        y[i] = z0(x[i])\n",
    "    #    else: y[i] = z1(x[i])\n",
    "\n",
    "    return y\n",
    "\n",
    "#let's try plotting the y values\n",
    "\n",
    "x_mm = np.linspace(0, 10, 50)\n",
    "y_mm = generate_dataset(-2, 4, x_mm)\n",
    "\n",
    "plt.scatter(x_mm, y_mm)\n",
    "plt.title('Synthetic Dataset')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d93bb271",
   "metadata": {},
   "source": [
    "### Problem 2b\n",
    "\n",
    "Write a function that implements the conditional density\n",
    "$$\n",
    "p(y | x; \\theta, \\beta)\n",
    "$$\n",
    "assuming that each $y^i$ is conditionally independent of each other given $z_i$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6dab3da4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixture_model_density(theta: float, beta: float, X: np.ndarray, y: np.ndarray) -> float:\n",
    "    # X is a one dimensional array of x values, N\n",
    "    # dim(y) = dim(x)\n",
    "\n",
    "    normals0 = []\n",
    "    normals1 = []\n",
    "\n",
    "    normals0 = [ stats.norm(loc=np.dot(theta, X[i]), scale=1) for i in range(X.shape[0]) ]\n",
    "    normals1 = [ stats.norm(loc=np.dot(beta, X[i]), scale=1) for i in range(X.shape[0]) ]\n",
    "\n",
    "    #for i in range(x.shape[0]):\n",
    "\n",
    "    return np.prod([(0.5 * normals0[i].pdf(y[i]) + 0.5 * normals1[i].pdf(y[i])) for i in range(len(y))])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6c57116",
   "metadata": {},
   "source": [
    "### Problem 2c\n",
    "\n",
    "Write a function that finds the optimal parameters of $\\theta$ and $\\beta$, i.e., \n",
    "$$\n",
    "\\operatorname{argmin}_{\\theta} p(y | x; \\theta, \\beta) \\,.\n",
    "$$\n",
    "You can use `generate_dataset` to test if your code is working."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "287fb2ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "63049e3b29464184a1bc614a4e0845e4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/50 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(theta, beta): (-1.8775510204081636, 4.0)\n",
      "Density based on these weights: 2.9332597770789448e-43\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(None, None)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def solve_mixture_model_density(X: np.ndarray, Y: np.ndarray) -> Tuple[float, float]:\n",
    "    # X contains a length N array of inputs of dimension 1\n",
    "    # Y contains a length N array of outputs\n",
    "    \n",
    "    best_theta = 0.0\n",
    "    best_beta = 0.0\n",
    "    best_density = 0.0\n",
    "        \n",
    "    searchspace = np.linspace(-4, 4)\n",
    "\n",
    "    for theta in tqdm(searchspace):\n",
    "        for beta in searchspace:\n",
    "            density = mixture_model_density(theta, beta, X, Y)\n",
    "            if density > best_density:\n",
    "                best_density = density\n",
    "                best_theta = theta\n",
    "                best_beta = beta\n",
    "    \n",
    "    return best_theta, best_beta\n",
    "\n",
    "thetas_mm = solve_mixture_model_density(x_mm, y_mm)\n",
    "\n",
    "print(\"(theta, beta): \" + str(thetas_mm)), print(\"Density based on these weights: \" + str(mixture_model_density(thetas_mm[0], thetas_mm[1], x_mm, y_mm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb03b443",
   "metadata": {},
   "source": [
    "## Problem 3: Linear Regression, Duplicate Data\n",
    "\n",
    "Suppose we have a dataset $(x_i, y_i)_{1 \\leq i \\leq N}$.\n",
    "Suppose we duplicate this dataset so that it now contains a copy of each item.\n",
    "Explain in words what happens to the weights of the linear regression in terms of the linear algebra solution and in terms of the probabilistic interpretation. Compare and contrast the two solutions. You can additionally write code and give experimental evidence as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00bb2cf1",
   "metadata": {},
   "source": [
    "### Answer: \n",
    "\n",
    "In terms of the linear algebra solution, the weights of the linear regression should not change when the dataset is duplicated. This is because the linear algebra solution is simply a matrix multiplication, and duplicate rows in the dataset do not change the product.\n",
    "\n",
    "In terms of the probabilistic interpretation, the weights of the linear regression should change when the dataset is duplicated. This is because the probabilistic interpretation is based on minimizing the sum of squared errors, and doubling the dataset doubles the sum of squared errors. Therefore, the weights that minimize the sum of squared errors will be different when the dataset is duplicated.\n",
    "\n",
    "Extra Notes: \n",
    "1) Weights stay the same in both cases.\n",
    "2) Loss increases with duplicate data.\n",
    "3) Mean (weights) stays the same, but the variance decreases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b56d3fa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function: 4x+5: using 1a, theta1 = 4, theta2 = 5\n",
    "\n",
    "x1 = [2., 1.]\n",
    "x2 = [4., 1.]\n",
    "x3 = [6., 1.]\n",
    "x4 = [8., 1.]\n",
    "x5 = [10., 1.]\n",
    "\n",
    "x = np.array([x1, x2, x3, x4, x5])\n",
    "y = [13., 21., 29., 37., 45.]\n",
    "\n",
    "# repeat the elements of x and y 5 times\n",
    "\n",
    "x5 = np.repeat(x, 5, axis=0)\n",
    "y5 = np.repeat(y, 5, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ebcb25dc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Thetas obtained from linear algebra (with duplications): [4. 5.]\n"
     ]
    }
   ],
   "source": [
    "thetas = np.linalg.solve(x5.transpose() @ x5, x5.transpose() @ y5)\n",
    "\n",
    "print(\"Thetas obtained from linear algebra (with duplications): \" + str(thetas))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "781bb9f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp=[]\n",
    "for i in range(len(coords)):\n",
    "    if i not in temp:\n",
    "        for j in range(i+1,len(coords)):\n",
    "            c1=coords[i]\n",
    "            c2=coords[j]\n",
    "            C1=SkyCoord(c1[0],c1[1], frame='galactic',unit=u.deg)\n",
    "            C2=SkyCoord(c2[0],c2[1], frame='galactic',unit=u.deg)\n",
    "            if C1.separation(C2)<=0.005*u.deg:\n",
    "                temp.append(j)\n",
    "new_coords=[]\n",
    "for i in range(len(coords)):\n",
    "    if i in temp:\n",
    "        continue\n",
    "    else:\n",
    "        new_coords.append(coords[i])\n",
    "coords=new_coords\n",
    "len(coords)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "40d3a090f54c6569ab1632332b64b2c03c39dcf918b08424e98f38b5ae0af88f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}