15. Множественная регрессия.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta [30.619881701311712, 0.9702056472470465, -1.8671913880379478, 0.9163711597955347]\n",
      "r-квадрат 0.6800074955952597\n",
      "\n",
      "отступление: бутстрапирование\n",
      "\n",
      "bootstrap_statistic(close_to_100, median, 100):\n",
      "[100.02634537429101, 100.06586584030461, 100.06850285070729, 100.06586584030461, 100.08022923777746, 100.02634537429101, 100.0251392489555, 100.06586584030461, 100.0251392489555, 100.18315165595742, 100.06627711829937, 100.06586584030461, 100.06627711829937, 100.06627711829937, 100.08022923777746, 100.08022923777746, 99.97928565593025, 99.95763666335068, 100.06627711829937, 99.95763666335068, 100.06586584030461, 100.06627711829937, 100.06850285070729, 100.08732646013036, 100.0340458696825, 100.02634537429101, 100.16675086620182, 100.11133261277095, 100.02634537429101, 100.0340458696825, 100.02048509159759, 100.06586584030461, 100.0340458696825, 100.02048509159759, 100.06850285070729, 100.06850285070729, 99.95763666335068, 100.08022923777746, 100.02048509159759, 100.07674830805912, 100.06586584030461, 100.06850285070729, 100.06627711829937, 100.11133261277095, 100.08022923777746, 100.06627711829937, 100.02634537429101, 100.0251392489555, 100.06850285070729, 100.06627711829937, 100.06586584030461, 100.0251392489555, 100.08022923777746, 100.11133261277095, 100.06586584030461, 100.0251392489555, 100.02634537429101, 100.04153611616675, 100.08022923777746, 100.06627711829937, 100.19440587880706, 100.04153611616675, 99.98693719710383, 100.06850285070729, 100.07674830805912, 100.08022923777746, 100.08732646013036, 100.06627711829937, 100.04153611616675, 100.08732646013036, 100.11133261277095, 100.06586584030461, 100.0251392489555, 100.13015471695907, 100.02634537429101, 99.95728292174014, 100.08732646013036, 100.06850285070729, 99.92424918869057, 100.02048509159759, 99.97928565593025, 99.92424918869057, 100.13015471695907, 100.06586584030461, 100.08022923777746, 100.06850285070729, 100.04153611616675, 100.06627711829937, 100.06586584030461, 100.06850285070729, 100.02048509159759, 100.02634537429101, 100.06627711829937, 100.04153611616675, 100.02634537429101, 100.02634537429101, 100.06627711829937, 100.18315165595742, 100.06850285070729, 100.06627711829937]\n",
      "\n",
      "bootstrap_statistic(far_from_100, median, 100):\n",
      "[0.9746929609488131, 0.9746929609488131, 0.9632764195933914, 200.0254523027627, 200.03723940374888, 0.9740127081653535, 200.0254523027627, 200.03723940374888, 0.9740127081653535, 200.043346407288, 200.06556682467556, 0.9679946120377706, 0.9229610720300252, 0.9170498145371947, 0.9746929609488131, 200.02972914757345, 100.37472916630605, 200.0439362129153, 0.9170498145371947, 0.9746929609488131, 0.9819299666825538, 0.9740127081653535, 200.02972914757345, 0.9819299666825538, 200.0439362129153, 200.03611456820246, 200.0439362129153, 200.02972914757345, 200.03611456820246, 200.0254523027627, 200.03723940374888, 0.9819299666825538, 0.9229610720300252, 0.9413891648644847, 0.9413891648644847, 200.01173496311017, 0.9746929609488131, 200.02972914757345, 200.043346407288, 0.9679946120377706, 200.07626904987512, 200.01173496311017, 0.9819299666825538, 200.01173496311017, 0.9632764195933914, 0.9746929609488131, 0.9746929609488131, 200.01173496311017, 0.9632764195933914, 0.9746929609488131, 200.02972914757345, 0.8082392572934448, 0.9413891648644847, 200.11842771777464, 0.9229610720300252, 0.9229610720300252, 100.37472916630605, 0.9819299666825538, 0.9740127081653535, 200.01173496311017, 200.03611456820246, 0.9819299666825538, 200.01173496311017, 0.9632764195933914, 200.0254523027627, 200.02972914757345, 200.02972914757345, 200.03611456820246, 0.9413891648644847, 100.37472916630605, 0.9170498145371947, 200.03723940374888, 0.9632764195933914, 0.8082392572934448, 0.9170498145371947, 200.01173496311017, 0.9746929609488131, 0.9746929609488131, 200.0439362129153, 100.37472916630605, 200.01173496311017, 200.03611456820246, 200.0254523027627, 0.9632764195933914, 200.02972914757345, 0.9740127081653535, 200.0254523027627, 100.37472916630605, 200.07626904987512, 0.9679946120377706, 200.02972914757345, 200.03723940374888, 200.01173496311017, 200.07626904987512, 0.9746929609488131, 0.9740127081653535, 0.9819299666825538, 200.01173496311017, 100.37472916630605, 100.37472916630605]\n",
      "\n",
      "\n",
      "стандартные ошибки бутстрапированных выборок [0.953551702104508, 0.06288763616183773, 0.11722269488203318, 0.8591786495949066]\n",
      "\n",
      "p_value(30.63, 1.174) 0.0\n",
      "p_value(0.972, 0.079) 0.0\n",
      "p_value(-1.868, 0.131) 0.0\n",
      "p_value(0.911, 0.990) 0.35746719881669264\n",
      "\n",
      "регуляризация\n",
      "\n",
      "alpha 0.0\n",
      "beta [30.619881701311712, 0.9702056472470465, -1.8671913880379478, 0.9163711597955347]\n",
      "скалярное_произведение(beta[1:],beta[1:]) 5.267438780018153\n",
      "r-квадрат 0.6800074955952597\n",
      "\n",
      "alpha 0.01\n",
      "beta [30.55985204967343, 0.9730655363505671, -1.8624424625144256, 0.9317665551046306]\n",
      "скалярное_произведение(beta[1:],beta[1:]) 5.2837373774215655\n",
      "r-квадрат 0.680010213297079\n",
      "\n",
      "alpha 0.1\n",
      "beta [30.894860179735474, 0.9490275238632391, -1.8501720889216575, 0.5325129720515789]\n",
      "скалярное_произведение(beta[1:],beta[1:]) 4.607360065077926\n",
      "r-квадрат 0.6797276241305292\n",
      "\n",
      "alpha 1\n",
      "beta [30.666778908554885, 0.908635996761392, -1.6938673046100265, 0.09370161190283018]\n",
      "скалярное_произведение(beta[1:],beta[1:]) 3.7035858123105934\n",
      "r-квадрат 0.6757061537631815\n",
      "\n",
      "alpha 10\n",
      "beta [28.372861060795607, 0.7307660860322116, -0.9212163182015426, -0.018495551723207087]\n",
      "скалярное_произведение(beta[1:],beta[1:]) 1.3830006628491893\n",
      "r-квадрат 0.5752138470466858\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "# multiple_regression.py\n",
    "\n",
    "import sys\n",
    "sys.path.append(\"../code-python3-ru\")\n",
    "\n",
    "from collections import Counter\n",
    "from functools import partial\n",
    "from lib.linear_algebra import dot, vector_add\n",
    "from lib.statistics import median, standard_deviation\n",
    "from lib.probability import normal_cdf\n",
    "from lib.gradient_descent import minimize_stochastic\n",
    "from lib.simple_linear_regression import total_sum_of_squares\n",
    "import math, random\n",
    "\n",
    "def predict(x_i, beta):\n",
    "    return dot(x_i, beta)\n",
    "\n",
    "def error(x_i, y_i, beta):\n",
    "    return y_i - predict(x_i, beta)\n",
    "\n",
    "def squared_error(x_i, y_i, beta):\n",
    "    return error(x_i, y_i, beta) ** 2\n",
    "\n",
    "def squared_error_gradient(x_i, y_i, beta):\n",
    "    \"\"\"градиент (относительно beta), соответствующий i-ой квадратичной ошибке\"\"\"\n",
    "    return [-2 * x_ij * error(x_i, y_i, beta)\n",
    "            for x_ij in x_i]\n",
    "\n",
    "def estimate_beta(x, y):\n",
    "    beta_initial = [random.random() for x_i in x[0]]\n",
    "    return minimize_stochastic(squared_error,\n",
    "                               squared_error_gradient,\n",
    "                               x, y,\n",
    "                               beta_initial,\n",
    "                               0.001)\n",
    "\n",
    "def multiple_r_squared(x, y, beta):\n",
    "    sum_of_squared_errors = sum(error(x_i, y_i, beta) ** 2\n",
    "                                for x_i, y_i in zip(x, y))\n",
    "    return 1.0 - sum_of_squared_errors / total_sum_of_squares(y)\n",
    "\n",
    "def bootstrap_sample(data):\n",
    "    \"\"\"случайно отбирает len(data) элементов с возвратом (с повторами\"\"\"\n",
    "    return [random.choice(data) for _ in data]\n",
    "\n",
    "def bootstrap_statistic(data, stats_fn, num_samples):\n",
    "    \"\"\"оценивает функцию stats_fn на бустрап-выборках из данных\n",
    "    в количестве num_samples\"\"\"\n",
    "    return [stats_fn(bootstrap_sample(data))\n",
    "            for _ in range(num_samples)]\n",
    "\n",
    "def estimate_sample_beta(sample):\n",
    "    x_sample, y_sample = list(zip(*sample)) # трюк с \"разъединением\" списка\n",
    "    return estimate_beta(x_sample, y_sample)\n",
    "\n",
    "def p_value(beta_hat_j, sigma_hat_j):\n",
    "    if beta_hat_j > 0:\n",
    "        return 2 * (1 - normal_cdf(beta_hat_j / sigma_hat_j))\n",
    "    else:\n",
    "        return 2 * normal_cdf(beta_hat_j / sigma_hat_j)\n",
    "\n",
    "#\n",
    "# РЕГУЛЯРИЗОВАННАЯ РЕГРЕССИЯ\n",
    "#\n",
    "\n",
    "# alpha - это ГИПЕРПАРАМЕТР, контролирующий, насколько жестким будет штраф\n",
    "# иногада его называют лямбда-параметром, но этот термин уже имеет свое\n",
    "# значение в языке Python\n",
    "def ridge_penalty(beta, alpha):\n",
    "  return alpha * dot(beta[1:], beta[1:])\n",
    "\n",
    "def squared_error_ridge(x_i, y_i, beta, alpha):\n",
    "    \"\"\"оценить ошибку плюс штраф для beta\"\"\"\n",
    "    return error(x_i, y_i, beta) ** 2 + ridge_penalty(beta, alpha)\n",
    "\n",
    "def ridge_penalty_gradient(beta, alpha):\n",
    "    \"\"\"градиент штрафа в гребневой регрессии\"\"\"\n",
    "    return [0] + [2 * alpha * beta_j for beta_j in beta[1:]]\n",
    "\n",
    "def squared_error_ridge_gradient(x_i, y_i, beta, alpha):\n",
    "    \"\"\"градиент, соответствующий i-ой квадратичной ошибке,\n",
    "    включая гребневой штраф\"\"\"\n",
    "    return vector_add(squared_error_gradient(x_i, y_i, beta),\n",
    "                      ridge_penalty_gradient(beta, alpha))\n",
    "\n",
    "def estimate_beta_ridge(x, y, alpha):\n",
    "    \"\"\"применить градиентный спуск для подгонки гребневой регрессии\n",
    "    со штрафом alpha\"\"\"\n",
    "    beta_initial = [random.random() for x_i in x[0]]\n",
    "    return minimize_stochastic(partial(squared_error_ridge, alpha=alpha),\n",
    "                               partial(squared_error_ridge_gradient,\n",
    "                                       alpha=alpha),\n",
    "                               x, y,\n",
    "                               beta_initial,\n",
    "                               0.001)\n",
    "\n",
    "def lasso_penalty(beta, alpha):\n",
    "    return alpha * sum(abs(beta_i) for beta_i in beta[1:])\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "\n",
    "    x = [[1,49,4,0],[1,41,9,0],[1,40,8,0],[1,25,6,0],[1,21,1,0],\n",
    "         [1,21,0,0],[1,19,3,0],[1,19,0,0],[1,18,9,0],[1,18,8,0],\n",
    "         [1,16,4,0],[1,15,3,0],[1,15,0,0],[1,15,2,0],[1,15,7,0],\n",
    "         [1,14,0,0],[1,14,1,0],[1,13,1,0],[1,13,7,0],[1,13,4,0],\n",
    "         [1,13,2,0],[1,12,5,0],[1,12,0,0],[1,11,9,0],[1,10,9,0],\n",
    "         [1,10,1,0],[1,10,1,0],[1,10,7,0],[1,10,9,0],[1,10,1,0],\n",
    "         [1,10,6,0],[1,10,6,0],[1,10,8,0],[1,10,10,0],[1,10,6,0],\n",
    "         [1,10,0,0],[1,10,5,0],[1,10,3,0],[1,10,4,0],[1,9,9,0],\n",
    "         [1,9,9,0],[1,9,0,0],[1,9,0,0],[1,9,6,0],[1,9,10,0],\n",
    "         [1,9,8,0],[1,9,5,0],[1,9,2,0],[1,9,9,0],[1,9,10,0],\n",
    "         [1,9,7,0],[1,9,2,0],[1,9,0,0],[1,9,4,0],[1,9,6,0],\n",
    "         [1,9,4,0],[1,9,7,0],[1,8,3,0],[1,8,2,0],[1,8,4,0],\n",
    "         [1,8,9,0],[1,8,2,0],[1,8,3,0],[1,8,5,0],[1,8,8,0],\n",
    "         [1,8,0,0],[1,8,9,0],[1,8,10,0],[1,8,5,0],[1,8,5,0],\n",
    "         [1,7,5,0],[1,7,5,0],[1,7,0,0],[1,7,2,0],[1,7,8,0],\n",
    "         [1,7,10,0],[1,7,5,0],[1,7,3,0],[1,7,3,0],[1,7,6,0],\n",
    "         [1,7,7,0],[1,7,7,0],[1,7,9,0],[1,7,3,0],[1,7,8,0],\n",
    "         [1,6,4,0],[1,6,6,0],[1,6,4,0],[1,6,9,0],[1,6,0,0],\n",
    "         [1,6,1,0],[1,6,4,0],[1,6,1,0],[1,6,0,0],[1,6,7,0],\n",
    "         [1,6,0,0],[1,6,8,0],[1,6,4,0],[1,6,2,1],[1,6,1,1],\n",
    "         [1,6,3,1],[1,6,6,1],[1,6,4,1],[1,6,4,1],[1,6,1,1],\n",
    "         [1,6,3,1],[1,6,4,1],[1,5,1,1],[1,5,9,1],[1,5,4,1],\n",
    "         [1,5,6,1],[1,5,4,1],[1,5,4,1],[1,5,10,1],[1,5,5,1],\n",
    "         [1,5,2,1],[1,5,4,1],[1,5,4,1],[1,5,9,1],[1,5,3,1],\n",
    "         [1,5,10,1],[1,5,2,1],[1,5,2,1],[1,5,9,1],[1,4,8,1],\n",
    "         [1,4,6,1],[1,4,0,1],[1,4,10,1],[1,4,5,1],[1,4,10,1],\n",
    "         [1,4,9,1],[1,4,1,1],[1,4,4,1],[1,4,4,1],[1,4,0,1],\n",
    "         [1,4,3,1],[1,4,1,1],[1,4,3,1],[1,4,2,1],[1,4,4,1],\n",
    "         [1,4,4,1],[1,4,8,1],[1,4,2,1],[1,4,4,1],[1,3,2,1],\n",
    "         [1,3,6,1],[1,3,4,1],[1,3,7,1],[1,3,4,1],[1,3,1,1],\n",
    "         [1,3,10,1],[1,3,3,1],[1,3,4,1],[1,3,7,1],[1,3,5,1],\n",
    "         [1,3,6,1],[1,3,1,1],[1,3,6,1],[1,3,10,1],[1,3,2,1],\n",
    "         [1,3,4,1],[1,3,2,1],[1,3,1,1],[1,3,5,1],[1,2,4,1],\n",
    "         [1,2,2,1],[1,2,8,1],[1,2,3,1],[1,2,1,1],[1,2,9,1],\n",
    "         [1,2,10,1],[1,2,9,1],[1,2,4,1],[1,2,5,1],[1,2,0,1],\n",
    "         [1,2,9,1],[1,2,9,1],[1,2,0,1],[1,2,1,1],[1,2,1,1],\n",
    "         [1,2,4,1],[1,1,0,1],[1,1,2,1],[1,1,2,1],[1,1,5,1],\n",
    "         [1,1,3,1],[1,1,10,1],[1,1,6,1],[1,1,0,1],[1,1,8,1],\n",
    "         [1,1,6,1],[1,1,4,1],[1,1,9,1],[1,1,9,1],[1,1,4,1],\n",
    "         [1,1,2,1],[1,1,9,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],\n",
    "         [1,1,1,1],[1,1,1,1],[1,1,5,1]]\n",
    "\n",
    "    daily_minutes_good = [68.77,51.25,52.08,38.36,44.54,57.13,51.4,41.42,31.22,34.76,\n",
    "                          54.01,38.79,47.59,49.1,27.66,41.03,36.73,48.65,28.12,46.62,\n",
    "                          35.57,32.98,35,26.07,23.77,39.73,40.57,31.65,31.21,36.32,\n",
    "                          20.45,21.93,26.02,27.34,23.49,46.94,30.5,33.8,24.23,21.4,\n",
    "                          27.94,32.24,40.57,25.07,19.42,22.39,18.42,46.96,23.72,26.41,\n",
    "                          26.97,36.76,40.32,35.02,29.47,30.2,31,38.11,38.18,36.31,21.03,\n",
    "                          30.86,36.07,28.66,29.08,37.28,15.28,24.17,22.31,30.17,25.53,\n",
    "                          19.85,35.37,44.6,17.23,13.47,26.33,35.02,32.09,24.81,19.33,\n",
    "                          28.77,24.26,31.98,25.73,24.86,16.28,34.51,15.23,39.72,40.8,\n",
    "                          26.06,35.76,34.76,16.13,44.04,18.03,19.65,32.62,35.59,39.43,\n",
    "                          14.18,35.24,40.13,41.82,35.45,36.07,43.67,24.61,20.9,21.9,\n",
    "                          18.79,27.61,27.21,26.61,29.77,20.59,27.53,13.82,33.2,25,33.1,\n",
    "                          36.65,18.63,14.87,22.2,36.81,25.53,24.62,26.25,18.21,28.08,\n",
    "                          19.42,29.79,32.8,35.99,28.32,27.79,35.88,29.06,36.28,14.1,\n",
    "                          36.63,37.49,26.9,18.58,38.48,24.48,18.95,33.55,14.24,29.04,\n",
    "                          32.51,25.63,22.22,19,32.73,15.16,13.9,27.2,32.01,29.27,33,\n",
    "                          13.74,20.42,27.32,18.23,35.35,28.48,9.08,24.62,20.12,35.26,\n",
    "                          19.92,31.02,16.49,12.16,30.7,31.22,34.65,13.13,27.51,33.2,\n",
    "                          31.57,14.1,33.42,17.44,10.12,24.42,9.82,23.39,30.93,15.03,\n",
    "                          21.67,31.09,33.29,22.61,26.89,23.48,8.38,27.81,32.35,23.84]\n",
    "    random.seed(0)\n",
    "    beta = estimate_beta(x, daily_minutes_good) # [30.63, 0.972, -1.868, 0.911]\n",
    "    print(\"beta\", beta)\n",
    "    print(\"r-квадрат\", multiple_r_squared(x, daily_minutes_good, beta))\n",
    "    print()\n",
    "\n",
    "    print(\"отступление: бутстрапирование\")\n",
    "    # 101 точек, чьи значения очень близки к 100\n",
    "    close_to_100 = [99.5 + random.random() for _ in range(101)]\n",
    "\n",
    "    # 101 точек, из них 50 значений близки к 0, а другие 50 близки к 200\n",
    "    far_from_100 = ([99.5 + random.random()] +\n",
    "                    [random.random() for _ in range(50)] +\n",
    "                    [200 + random.random() for _ in range(50)])\n",
    "\n",
    "    print()\n",
    "    print(\"bootstrap_statistic(close_to_100, median, 100):\")\n",
    "    print(bootstrap_statistic(close_to_100, median, 100))\n",
    "    print()\n",
    "    print(\"bootstrap_statistic(far_from_100, median, 100):\")\n",
    "    print(bootstrap_statistic(far_from_100, median, 100))\n",
    "    print()\n",
    "\n",
    "    random.seed(0) # чтобы получить повторимые результаты\n",
    "\n",
    "    bootstrap_betas = bootstrap_statistic(list(zip(x, daily_minutes_good)),\n",
    "                                          estimate_sample_beta,\n",
    "                                          100)\n",
    "\n",
    "    bootstrap_standard_errors = [\n",
    "        standard_deviation([beta[i] for beta in bootstrap_betas])\n",
    "        for i in range(4)]\n",
    "\n",
    "    print()\n",
    "    print(\"стандартные ошибки бутстрапированных выборок\", bootstrap_standard_errors)\n",
    "    print()\n",
    "\n",
    "    print(\"p_value(30.63, 1.174)\", p_value(30.63, 1.174))\n",
    "    print(\"p_value(0.972, 0.079)\", p_value(0.972, 0.079))\n",
    "    print(\"p_value(-1.868, 0.131)\", p_value(-1.868, 0.131))\n",
    "    print(\"p_value(0.911, 0.990)\", p_value(0.911, 0.990))\n",
    "    print()\n",
    "\n",
    "    print(\"регуляризация\")\n",
    "    print()\n",
    "    \n",
    "    random.seed(0)\n",
    "    for alpha in [0.0, 0.01, 0.1, 1, 10]:\n",
    "        beta = estimate_beta_ridge(x, daily_minutes_good, alpha=alpha)\n",
    "        print(\"alpha\", alpha)\n",
    "        print(\"beta\", beta)\n",
    "        print(\"скалярное_произведение(beta[1:],beta[1:])\", dot(beta[1:], beta[1:]))\n",
    "        print(\"r-квадрат\", multiple_r_squared(x, daily_minutes_good, beta))\n",
    "        print()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}