Linear Fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "678762ac",
   "metadata": {},
   "source": [
    "# Linear Fitting\n",
    "\n",
    "**Problem:** Find `a` and `b` in the equation `y=a+b×x`, given the vectors `x` and `y`, using the [least-squares method](https://en.wikipedia.org/wiki/Linear_least_squares). Also, determine the corresponding [R-squared](https://en.wikipedia.org/wiki/Coefficient_of_determination) value."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "565a06e3",
   "metadata": {},
   "source": [
    "Test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b65828ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<svg version=\"1.1\" baseProfile=\"full\" width=\"100%\" height=\"100%\" viewBox=\"0 0 4320 3240\"  xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" >\n",
       "<desc>Created by Causeway SVG engine - SharpPlot v3.71.0</desc>\n",
       "<g stroke-miterlimit=\"100\">\n",
       "<desc>Paint the paper =====</desc>\n",
       " <g fill=\"White\" >\n",
       "  <rect x=\"0\" y=\"0\" width=\"4320\" height=\"3240\" />\n",
       " </g>\n",
       "<desc>Border =====</desc>\n",
       "<desc>Region =====</desc>\n",
       "<desc>X-Axis Ticks =====</desc>\n",
       "<desc>X-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L400 2764M685 2720L685 2764M969 2720L969 2764M1254 2720L1254 2764M1538 2720L1538 2764M1823 2720L1823 2764M2108 2720L2108 2764M2392 2720L2392 2764M2677 2720L2677 2764M2962 2720L2962 2764M3246 2720L3246 2764M3531 2720L3531 2764M3815 2720L3815 2764M4100 2720L4100 2764\"/>\n",
       "<desc>Y-Axis Ticks =====</desc>\n",
       "<desc>Y-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L367 2720M400 2546L367 2546M400 2372L367 2372M400 2198L367 2198M400 2025L367 2025M400 1851L367 1851M400 1677L367 1677M400 1503L367 1503M400 1329L367 1329M400 1155L367 1155M400 982L367 982M400 808L367 808M400 634L367 634M400 460L367 460\"/>\n",
       "<desc>Axes =====</desc>\n",
       " <path fill=\"none\" stroke-width=\"12\" stroke=\"Black\" d=\"M400 2764L400 416M356 2720L4144 2720\"/>\n",
       "<desc>Y-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"end\" >\n",
       "  <text x=\"335\" y=\"2740\" >10</text>\n",
       "  <text x=\"335\" y=\"2566\" >20</text>\n",
       "  <text x=\"335\" y=\"2392\" >30</text>\n",
       "  <text x=\"335\" y=\"2218\" >40</text>\n",
       "  <text x=\"335\" y=\"2045\" >50</text>\n",
       "  <text x=\"335\" y=\"1871\" >60</text>\n",
       "  <text x=\"335\" y=\"1697\" >70</text>\n",
       "  <text x=\"335\" y=\"1523\" >80</text>\n",
       "  <text x=\"335\" y=\"1349\" >90</text>\n",
       "  <text x=\"335\" y=\"1175\" >100</text>\n",
       "  <text x=\"335\" y=\"1002\" >110</text>\n",
       "  <text x=\"335\" y=\"828\" >120</text>\n",
       "  <text x=\"335\" y=\"654\" >130</text>\n",
       "  <text x=\"335\" y=\"480\" >140</text>\n",
       " </g>\n",
       "<desc>for X-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"middle\" >\n",
       "  <text x=\"400\" y=\"2844\" >0</text>\n",
       "  <text x=\"685\" y=\"2844\" >2</text>\n",
       "  <text x=\"969\" y=\"2844\" >4</text>\n",
       "  <text x=\"1254\" y=\"2844\" >6</text>\n",
       "  <text x=\"1538\" y=\"2844\" >8</text>\n",
       "  <text x=\"1823\" y=\"2844\" >10</text>\n",
       "  <text x=\"2108\" y=\"2844\" >12</text>\n",
       "  <text x=\"2392\" y=\"2844\" >14</text>\n",
       "  <text x=\"2677\" y=\"2844\" >16</text>\n",
       "  <text x=\"2962\" y=\"2844\" >18</text>\n",
       "  <text x=\"3246\" y=\"2844\" >20</text>\n",
       "  <text x=\"3531\" y=\"2844\" >22</text>\n",
       "  <text x=\"3815\" y=\"2844\" >24</text>\n",
       "  <text x=\"4100\" y=\"2844\" >26</text>\n",
       " </g>\n",
       "<desc>Heading, subheading and footnotes =====</desc>\n",
       "<desc>Start of Line Chart ===========</desc>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Navy\" stroke-linejoin=\"round\" stroke-linecap=\"round\" fill=\"none\" >\n",
       "  <path d=\"M3958 547L3815 634 3673 721 3531 808 3388 895 3246 982 3104 1068 2962 1155 2819 1242 2677 1329 2535 1416 2392 1503 2250 1590 2108 1677 1965 1764 1823 1851 1681 1938 1538 2025 1396 2112 1254 2198 1112 2285 969 2372 827 2459 685 2546 542 2633\" />\n",
       " </g>\n",
       "<desc>Reset to original origin</desc>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x←⍳25\n",
    "y←10+5×x\n",
    "]plot y x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6a568dd",
   "metadata": {},
   "source": [
    "This data is too perfect. Let's add some noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e0e18071",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<svg version=\"1.1\" baseProfile=\"full\" width=\"100%\" height=\"100%\" viewBox=\"0 0 4320 3240\"  xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" >\n",
       "<desc>Created by Causeway SVG engine - SharpPlot v3.71.0</desc>\n",
       "<g stroke-miterlimit=\"100\">\n",
       "<desc>Paint the paper =====</desc>\n",
       " <g fill=\"White\" >\n",
       "  <rect x=\"0\" y=\"0\" width=\"4320\" height=\"3240\" />\n",
       " </g>\n",
       "<desc>Border =====</desc>\n",
       "<desc>Region =====</desc>\n",
       "<desc>X-Axis Ticks =====</desc>\n",
       "<desc>X-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L400 2764M685 2720L685 2764M969 2720L969 2764M1254 2720L1254 2764M1538 2720L1538 2764M1823 2720L1823 2764M2108 2720L2108 2764M2392 2720L2392 2764M2677 2720L2677 2764M2962 2720L2962 2764M3246 2720L3246 2764M3531 2720L3531 2764M3815 2720L3815 2764M4100 2720L4100 2764\"/>\n",
       "<desc>Y-Axis Ticks =====</desc>\n",
       "<desc>Y-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L367 2720M400 2546L367 2546M400 2372L367 2372M400 2198L367 2198M400 2025L367 2025M400 1851L367 1851M400 1677L367 1677M400 1503L367 1503M400 1329L367 1329M400 1155L367 1155M400 982L367 982M400 808L367 808M400 634L367 634M400 460L367 460\"/>\n",
       "<desc>Axes =====</desc>\n",
       " <path fill=\"none\" stroke-width=\"12\" stroke=\"Black\" d=\"M400 2764L400 416M356 2720L4144 2720\"/>\n",
       "<desc>Y-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"end\" >\n",
       "  <text x=\"335\" y=\"2740\" >10</text>\n",
       "  <text x=\"335\" y=\"2566\" >20</text>\n",
       "  <text x=\"335\" y=\"2392\" >30</text>\n",
       "  <text x=\"335\" y=\"2218\" >40</text>\n",
       "  <text x=\"335\" y=\"2045\" >50</text>\n",
       "  <text x=\"335\" y=\"1871\" >60</text>\n",
       "  <text x=\"335\" y=\"1697\" >70</text>\n",
       "  <text x=\"335\" y=\"1523\" >80</text>\n",
       "  <text x=\"335\" y=\"1349\" >90</text>\n",
       "  <text x=\"335\" y=\"1175\" >100</text>\n",
       "  <text x=\"335\" y=\"1002\" >110</text>\n",
       "  <text x=\"335\" y=\"828\" >120</text>\n",
       "  <text x=\"335\" y=\"654\" >130</text>\n",
       "  <text x=\"335\" y=\"480\" >140</text>\n",
       " </g>\n",
       "<desc>for X-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"middle\" >\n",
       "  <text x=\"400\" y=\"2844\" >0</text>\n",
       "  <text x=\"685\" y=\"2844\" >2</text>\n",
       "  <text x=\"969\" y=\"2844\" >4</text>\n",
       "  <text x=\"1254\" y=\"2844\" >6</text>\n",
       "  <text x=\"1538\" y=\"2844\" >8</text>\n",
       "  <text x=\"1823\" y=\"2844\" >10</text>\n",
       "  <text x=\"2108\" y=\"2844\" >12</text>\n",
       "  <text x=\"2392\" y=\"2844\" >14</text>\n",
       "  <text x=\"2677\" y=\"2844\" >16</text>\n",
       "  <text x=\"2962\" y=\"2844\" >18</text>\n",
       "  <text x=\"3246\" y=\"2844\" >20</text>\n",
       "  <text x=\"3531\" y=\"2844\" >22</text>\n",
       "  <text x=\"3815\" y=\"2844\" >24</text>\n",
       "  <text x=\"4100\" y=\"2844\" >26</text>\n",
       " </g>\n",
       "<desc>Heading, subheading and footnotes =====</desc>\n",
       "<desc>Start of Line Chart ===========</desc>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Navy\" stroke-linejoin=\"round\" stroke-linecap=\"round\" fill=\"none\" >\n",
       "  <path d=\"M3958 610L3815 691 3673 769 3531 820 3388 974 3246 1011 3104 1104 2962 1159 2819 1281 2677 1331 2535 1395 2392 1521 2250 1596 2108 1714 1965 1708 1823 1927 1681 1914 1538 2086 1396 2140 1254 2191 1112 2307 969 2404 827 2526 685 2541 542 2698\" />\n",
       " </g>\n",
       "<desc>Reset to original origin</desc>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y+←(5-⍨10×?0⍴⍨≢y)\n",
    "]plot y x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa0151da",
   "metadata": {},
   "source": [
    "The slope is easy to determine if we already know what is the value of `a`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b990853f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">4.900895018\n",
       "</pre>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "⊢b←(y-10)⌹x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dfd85f5",
   "metadata": {},
   "source": [
    "We can also use `⌹` to fit several variables. To fit the value of `a`, we consider an additional variable with a constant value of 1, as explained in the [docuentation for `⌹`](https://help.dyalog.com/latest/index.htm#Language/Primitive%20Functions/Matrix%20Divide.htm). We just need to modify our right argument from"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "98eca846",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25\n",
       "</pre>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d529781c",
   "metadata": {},
   "source": [
    "into"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c94d2abd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">1  1\n",
       "1  2\n",
       "1  3\n",
       "1  4\n",
       "1  5\n",
       "1  6\n",
       "1  7\n",
       "1  8\n",
       "1  9\n",
       "1 10\n",
       "1 11\n",
       "1 12\n",
       "1 13\n",
       "1 14\n",
       "1 15\n",
       "1 16\n",
       "1 17\n",
       "1 18\n",
       "1 19\n",
       "1 20\n",
       "1 21\n",
       "1 22\n",
       "1 23\n",
       "1 24\n",
       "1 25\n",
       "</pre>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1,⍪x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a99101d6",
   "metadata": {},
   "source": [
    "(a vector is treated by `⌹` as a column vector). Then:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5dbfe77b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">8.968943213 4.961545417\n",
       "</pre>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "⊢(a b)←y⌹1,⍪x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b2a9d64",
   "metadata": {},
   "source": [
    "## Fitting function\n",
    "\n",
    "Let's write a function to perform linear fittings. Our function will take as right argument the `y x` data and will be ambivalent (we use `0=⎕NC'⍺'` to check if a left argument was passed to the function). If there is a left argument, it will be interpreted as the value of the ordinate at the origin. So, we will fit the slope of the line through the origin using the data `⍵-⍺ 0`. If no argument is given, `a` and `b` will be fitted. In both cases, the function will return the values `a b`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d230f9bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">10 4.900895018\n",
       "</pre>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">8.968943213 4.961545417\n",
       "</pre>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<svg version=\"1.1\" baseProfile=\"full\" width=\"100%\" height=\"100%\" viewBox=\"0 0 4320 3240\"  xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" >\n",
       "<desc>Created by Causeway SVG engine - SharpPlot v3.71.0</desc>\n",
       "<g stroke-miterlimit=\"100\">\n",
       "<desc>Paint the paper =====</desc>\n",
       " <g fill=\"White\" >\n",
       "  <rect x=\"0\" y=\"0\" width=\"4320\" height=\"3240\" />\n",
       " </g>\n",
       "<desc>Border =====</desc>\n",
       "<desc>Region =====</desc>\n",
       "<desc>X-Axis Ticks =====</desc>\n",
       "<desc>X-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L400 2764M685 2720L685 2764M969 2720L969 2764M1254 2720L1254 2764M1538 2720L1538 2764M1823 2720L1823 2764M2108 2720L2108 2764M2392 2720L2392 2764M2677 2720L2677 2764M2962 2720L2962 2764M3246 2720L3246 2764M3531 2720L3531 2764M3815 2720L3815 2764M4100 2720L4100 2764\"/>\n",
       "<desc>Y-Axis Ticks =====</desc>\n",
       "<desc>Y-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L367 2720M400 2546L367 2546M400 2372L367 2372M400 2198L367 2198M400 2025L367 2025M400 1851L367 1851M400 1677L367 1677M400 1503L367 1503M400 1329L367 1329M400 1155L367 1155M400 982L367 982M400 808L367 808M400 634L367 634M400 460L367 460\"/>\n",
       "<desc>Axes =====</desc>\n",
       " <path fill=\"none\" stroke-width=\"12\" stroke=\"Black\" d=\"M400 2764L400 416M356 2720L4144 2720\"/>\n",
       "<desc>Y-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"end\" >\n",
       "  <text x=\"335\" y=\"2740\" >10</text>\n",
       "  <text x=\"335\" y=\"2566\" >20</text>\n",
       "  <text x=\"335\" y=\"2392\" >30</text>\n",
       "  <text x=\"335\" y=\"2218\" >40</text>\n",
       "  <text x=\"335\" y=\"2045\" >50</text>\n",
       "  <text x=\"335\" y=\"1871\" >60</text>\n",
       "  <text x=\"335\" y=\"1697\" >70</text>\n",
       "  <text x=\"335\" y=\"1523\" >80</text>\n",
       "  <text x=\"335\" y=\"1349\" >90</text>\n",
       "  <text x=\"335\" y=\"1175\" >100</text>\n",
       "  <text x=\"335\" y=\"1002\" >110</text>\n",
       "  <text x=\"335\" y=\"828\" >120</text>\n",
       "  <text x=\"335\" y=\"654\" >130</text>\n",
       "  <text x=\"335\" y=\"480\" >140</text>\n",
       " </g>\n",
       "<desc>for X-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"middle\" >\n",
       "  <text x=\"400\" y=\"2844\" >0</text>\n",
       "  <text x=\"685\" y=\"2844\" >2</text>\n",
       "  <text x=\"969\" y=\"2844\" >4</text>\n",
       "  <text x=\"1254\" y=\"2844\" >6</text>\n",
       "  <text x=\"1538\" y=\"2844\" >8</text>\n",
       "  <text x=\"1823\" y=\"2844\" >10</text>\n",
       "  <text x=\"2108\" y=\"2844\" >12</text>\n",
       "  <text x=\"2392\" y=\"2844\" >14</text>\n",
       "  <text x=\"2677\" y=\"2844\" >16</text>\n",
       "  <text x=\"2962\" y=\"2844\" >18</text>\n",
       "  <text x=\"3246\" y=\"2844\" >20</text>\n",
       "  <text x=\"3531\" y=\"2844\" >22</text>\n",
       "  <text x=\"3815\" y=\"2844\" >24</text>\n",
       "  <text x=\"4100\" y=\"2844\" >26</text>\n",
       " </g>\n",
       "<desc>Heading, subheading and footnotes =====</desc>\n",
       "<desc>Start of Line Chart ===========</desc>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Navy\" stroke-linejoin=\"round\" stroke-linecap=\"round\" fill=\"none\" >\n",
       "  <path d=\"M3958 590L3815 675 3673 760 3531 846 3388 931 3246 1016 3104 1101 2962 1186 2819 1272 2677 1357 2535 1442 2392 1527 2250 1612 2108 1698 1965 1783 1823 1868 1681 1953 1538 2038 1396 2124 1254 2209 1112 2294 969 2379 827 2464 685 2550 542 2635\" />\n",
       " </g>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Maroon\" stroke-linejoin=\"round\" stroke-linecap=\"round\" stroke-dasharray=\"9 12\" fill=\"none\" >\n",
       "  <path d=\"M3958 582L3815 668 3673 754 3531 840 3388 927 3246 1013 3104 1099 2962 1185 2819 1272 2677 1358 2535 1444 2392 1530 2250 1617 2108 1703 1965 1789 1823 1875 1681 1962 1538 2048 1396 2134 1254 2220 1112 2307 969 2393 827 2479 685 2565 542 2652\" />\n",
       " </g>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Teal\" stroke-linejoin=\"round\" stroke-linecap=\"round\" stroke-dasharray=\"0.3 5.7\" fill=\"none\" >\n",
       "  <path d=\"M3958 610L3815 691 3673 769 3531 820 3388 974 3246 1011 3104 1104 2962 1159 2819 1281 2677 1331 2535 1395 2392 1521 2250 1596 2108 1714 1965 1708 1823 1927 1681 1914 1538 2086 1396 2140 1254 2191 1112 2307 969 2404 827 2526 685 2541 542 2698\" />\n",
       " </g>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Green\" stroke-linejoin=\"round\" stroke-linecap=\"round\" stroke-dasharray=\"9 6.45 0.3 6.45\" fill=\"none\" >\n",
       "  <path d=\"M3958 610L3815 691 3673 769 3531 820 3388 974 3246 1011 3104 1104 2962 1159 2819 1281 2677 1331 2535 1395 2392 1521 2250 1596 2108 1714 1965 1708 1823 1927 1681 1914 1538 2086 1396 2140 1254 2191 1112 2307 969 2404 827 2526 685 2541 542 2698\" />\n",
       " </g>\n",
       "<desc>Reset to original origin</desc>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LF←{0=⎕NC'⍺': ⊃⌹∘(1,⍪)/⍵ ⋄ ⍺,⌹/⍵-⍺ 0}   ⍝ linear fit. ⍵: y x data. ⍺ (optional): ordinate\n",
    "⍝ eg\n",
    "⊢(ad bd)←10 LF y x\n",
    "⊢(a  b )←   LF y x\n",
    "]plot ((ad+bd×x) (a+b×x) y y) x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ac77307",
   "metadata": {},
   "source": [
    "## R-squared value\n",
    "\n",
    "We still need to determine the [R-squared](https://en.wikipedia.org/wiki/Coefficient_of_determination) value, which can be calculated as one minus the ratio of the sums of squares (`÷⍥(+/×⍨)`) of differences between the real and estimated values and between the real values and their mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0b3b5046",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">0.9971435494\n",
       "</pre>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1-(y-a+b×x)÷⍥(+/×⍨)(⊢-+⌿÷≢)y   ⍝ R-squared value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11f84fd",
   "metadata": {},
   "source": [
    "The R-squared value is defined not only for linear fittings, but also for any other fitting we perform on some data. Therefore, we can define an operator that takes the fitted function as operand (in our case, `{a+b×⍵}` or `a+b∘×`) and the fitting data as arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bf26cbd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">0.9971435494\n",
       "</pre>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_R2←{1-(⍺-⍺⍺⍵)÷⍥(+/×⍨)(⊢-+/÷≢)⍺}   ⍝ returns R-squared value for function ⍺⍺. ⍺ ⍵ is the y x data\n",
    "⍝ eg\n",
    "y (a+b∘×)_R2 x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4243e96a",
   "metadata": {},
   "source": [
    "More generally, we can define a dyadic operator `_R2_` that directly takes the fitting and fitted functions and returns the fitted parameters and R-squared value for the given data (as right argument) and optional parameters for the fitting function (as left argument). We then use this operator to define an `LR2` function that performs a linear fitting and returns the fitted parameters and the R-squared value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c895bec6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">10 4.900895018 0.996948701\n",
       "</pre>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">0 5.489130312 0.9823995937\n",
       "</pre>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<pre class=\"language-APL\">8.968943213 4.961545417 0.9971435494\n",
       "</pre>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<svg version=\"1.1\" baseProfile=\"full\" width=\"100%\" height=\"100%\" viewBox=\"0 0 4320 3240\"  xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" >\n",
       "<desc>Created by Causeway SVG engine - SharpPlot v3.71.0</desc>\n",
       "<g stroke-miterlimit=\"100\">\n",
       "<desc>Paint the paper =====</desc>\n",
       " <g fill=\"White\" >\n",
       "  <rect x=\"0\" y=\"0\" width=\"4320\" height=\"3240\" />\n",
       " </g>\n",
       "<desc>Border =====</desc>\n",
       "<desc>Region =====</desc>\n",
       "<desc>X-Axis Ticks =====</desc>\n",
       "<desc>X-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L400 2764M685 2720L685 2764M969 2720L969 2764M1254 2720L1254 2764M1538 2720L1538 2764M1823 2720L1823 2764M2108 2720L2108 2764M2392 2720L2392 2764M2677 2720L2677 2764M2962 2720L2962 2764M3246 2720L3246 2764M3531 2720L3531 2764M3815 2720L3815 2764M4100 2720L4100 2764\"/>\n",
       "<desc>Y-Axis Ticks =====</desc>\n",
       "<desc>Y-Axis tickmarks</desc>\n",
       " <path fill=\"none\" stroke-width=\"7.2\" stroke=\"Black\" d=\"M400 2720L367 2720M400 2559L367 2559M400 2397L367 2397M400 2236L367 2236M400 2074L367 2074M400 1913L367 1913M400 1751L367 1751M400 1590L367 1590M400 1429L367 1429M400 1267L367 1267M400 1106L367 1106M400 944L367 944M400 783L367 783M400 621L367 621M400 460L367 460\"/>\n",
       "<desc>Axes =====</desc>\n",
       " <path fill=\"none\" stroke-width=\"12\" stroke=\"Black\" d=\"M400 2764L400 416M356 2720L4144 2720\"/>\n",
       "<desc>Y-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"end\" >\n",
       "  <text x=\"335\" y=\"2740\" >0</text>\n",
       "  <text x=\"335\" y=\"2579\" >10</text>\n",
       "  <text x=\"335\" y=\"2417\" >20</text>\n",
       "  <text x=\"335\" y=\"2256\" >30</text>\n",
       "  <text x=\"335\" y=\"2094\" >40</text>\n",
       "  <text x=\"335\" y=\"1933\" >50</text>\n",
       "  <text x=\"335\" y=\"1771\" >60</text>\n",
       "  <text x=\"335\" y=\"1610\" >70</text>\n",
       "  <text x=\"335\" y=\"1449\" >80</text>\n",
       "  <text x=\"335\" y=\"1287\" >90</text>\n",
       "  <text x=\"335\" y=\"1126\" >100</text>\n",
       "  <text x=\"335\" y=\"964\" >110</text>\n",
       "  <text x=\"335\" y=\"803\" >120</text>\n",
       "  <text x=\"335\" y=\"641\" >130</text>\n",
       "  <text x=\"335\" y=\"480\" >140</text>\n",
       " </g>\n",
       "<desc>for X-axis labels</desc>\n",
       " <g font-family=\"Times New Roman\" font-size=\"80\" text-anchor=\"middle\" >\n",
       "  <text x=\"400\" y=\"2844\" >0</text>\n",
       "  <text x=\"685\" y=\"2844\" >2</text>\n",
       "  <text x=\"969\" y=\"2844\" >4</text>\n",
       "  <text x=\"1254\" y=\"2844\" >6</text>\n",
       "  <text x=\"1538\" y=\"2844\" >8</text>\n",
       "  <text x=\"1823\" y=\"2844\" >10</text>\n",
       "  <text x=\"2108\" y=\"2844\" >12</text>\n",
       "  <text x=\"2392\" y=\"2844\" >14</text>\n",
       "  <text x=\"2677\" y=\"2844\" >16</text>\n",
       "  <text x=\"2962\" y=\"2844\" >18</text>\n",
       "  <text x=\"3246\" y=\"2844\" >20</text>\n",
       "  <text x=\"3531\" y=\"2844\" >22</text>\n",
       "  <text x=\"3815\" y=\"2844\" >24</text>\n",
       "  <text x=\"4100\" y=\"2844\" >26</text>\n",
       " </g>\n",
       "<desc>Heading, subheading and footnotes =====</desc>\n",
       "<desc>Start of Line Chart ===========</desc>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Navy\" stroke-linejoin=\"round\" stroke-linecap=\"round\" fill=\"none\" >\n",
       "  <path d=\"M3958 581L3815 660 3673 739 3531 818 3388 897 3246 976 3104 1055 2962 1135 2819 1214 2677 1293 2535 1372 2392 1451 2250 1530 2108 1609 1965 1688 1823 1767 1681 1847 1538 1926 1396 2005 1254 2084 1112 2163 969 2242 827 2321 685 2400 542 2479\" />\n",
       " </g>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Maroon\" stroke-linejoin=\"round\" stroke-linecap=\"round\" stroke-dasharray=\"9 12\" fill=\"none\" >\n",
       "  <path d=\"M3958 505L3815 593 3673 682 3531 771 3388 859 3246 948 3104 1036 2962 1125 2819 1214 2677 1302 2535 1391 2392 1479 2250 1568 2108 1657 1965 1745 1823 1834 1681 1923 1538 2011 1396 2100 1254 2188 1112 2277 969 2366 827 2454 685 2543 542 2631\" />\n",
       " </g>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Teal\" stroke-linejoin=\"round\" stroke-linecap=\"round\" stroke-dasharray=\"0.3 5.7\" fill=\"none\" >\n",
       "  <path d=\"M3958 573L3815 653 3673 733 3531 813 3388 893 3246 973 3104 1053 2962 1134 2819 1214 2677 1294 2535 1374 2392 1454 2250 1534 2108 1614 1965 1694 1823 1774 1681 1854 1538 1934 1396 2015 1254 2095 1112 2175 969 2255 827 2335 685 2415 542 2495\" />\n",
       " </g>\n",
       "<desc>Points follow ...</desc>\n",
       "<desc>Line</desc>\n",
       " <g stroke-width=\"3\" stroke=\"Green\" stroke-linejoin=\"round\" stroke-linecap=\"round\" stroke-dasharray=\"9 6.45 0.3 6.45\" fill=\"none\" >\n",
       "  <path d=\"M3958 599L3815 675 3673 747 3531 795 3388 937 3246 972 3104 1058 2962 1109 2819 1222 2677 1269 2535 1328 2392 1445 2250 1515 2108 1624 1965 1619 1823 1822 1681 1810 1538 1970 1396 2020 1254 2067 1112 2175 969 2265 827 2378 685 2392 542 2538\" />\n",
       " </g>\n",
       "<desc>Reset to original origin</desc>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_R2_←{⍺←⊢ ⋄ r←⍺ ⍺⍺ ⍵ ⋄ r,r∘⍵⍵_R2/⍵}   ⍝ returns fitted parameters and R-squared value\n",
    "                                      ⍝ ⍺⍺: fitting function, ⍵⍵: fitted function\n",
    "                                      ⍝ ⍵: fitting data, ⍺ (optional): left argument of fitting function\n",
    "LR2←LF _R2_{(a b)←⍺ ⋄ (a+b∘×)⍵}       ⍝ linear fitting with y x data ⍵ and (optional) ordinate ⍺\n",
    "⍝ eg\n",
    "⊢ad bd _←10 LR2 y x\n",
    "⊢az bz _← 0 LR2 y x\n",
    "⊢a  b  _←   LR2 y x\n",
    "]plot ((ad+bd×x) (az+bz×x) (a+b×x) y) x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18dcc0c4",
   "metadata": {},
   "source": [
    "Obviously, the R-squared value is higher when we fit both values, since the fitted function can accomodate better the fitted data. However, it is useful to have an optional argument to impose already known values, as it will often be the case for functions that are known to pass by the origin."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0331fa15",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fd39b740-05e8-4534-b493-ffd0089fbf1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "_R2←{1-(⍺-⍺⍺⍵)÷⍥(+/×⍨)(⊢-+/÷≢)⍺}        ⍝ R-squared value for ⍺⍺ fitted function\n",
    "_R2_←{⍺←⊢ ⋄ r←⍺ ⍺⍺ ⍵ ⋄ r,r∘⍵⍵_R2/⍵}     ⍝ fit and R-squared value for fitting and fitted functions ⍺⍺ and ⍵⍵\n",
    "LF←{0=⎕NC'⍺': ⊃⌹∘(1,⍪)/⍵ ⋄ ⍺,⌹/⍵-⍺ 0}   ⍝ linear fitting function (⍵: y x data, ⍺: optional ordinate)\n",
    "LR2←LF _R2_{(a b)←⍺ ⋄ (a+b∘×)⍵}         ⍝ linear fitting and R-squared value (⍵: y x data, ⍺: optional ordinate)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Dyalog APL",
   "language": "apl",
   "name": "dyalog-kernel"
  },
  "language_info": {
   "file_extension": ".apl",
   "mimetype": "text/apl",
   "name": "APL"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}