{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A very simple example of fitting a function or model to data by using least squares.\n",
    "\n",
    "#The Data: Here are some hypothetical data in the form of (t,y) pairs:\n",
    "data = [[1.1, 1.24], [1.9, 0.83], [2.3, 0.71], [4.1, 0.29], [5.5, 0.15]];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot:\n",
    "plt1 = scatter_plot(data);\n",
    "show(plt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Model and Sum of Squares: Let's fit a model f(t) = a*exp(b*t) to this\n",
    "#data by adjusting a and b. First form a sum of squares SS(a,b):\n",
    "var('t a b');\n",
    "f(t,a,b) = a*exp(b*t); #The model function to fit\n",
    "SS = function('SS')(a,b);\n",
    "SS(a,b) = add((f(data[i][0],a,b)-data[i][1])^2 for i in range(5));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#With only two parameters \"a\" and \"b\", a visual estimate of the best choice\n",
    "#(the choices of a and b that minimize SS) can be found by plotting. It's\n",
    "#clear that b<0 since the data decays, and also that a>1. Let us plot log(SS):\n",
    "#A quick plot of log(SS(a,b))\n",
    "plot3d(log(SS), (a,0,3), (b,-1,0)).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Rotating the graph around shows a minimum around a = 2, b = -0.5.\n",
    "#We could minimize using Calculus 3 techniques, but Sage has a\n",
    "#built-in command \"find_fit\" to find the least-squares model fit to the data:\n",
    "model(t) = f(t,a,b); #Specify the model to be fit, with \"t\" as the independent variable\n",
    "sol = find_fit(data,model,parameters=[a,b],initial_guess=[2.0,-0.5]) #Fit the model by adjusting k and b\n",
    "abest = sol[0].rhs();\n",
    "bbest = sol[1].rhs();\n",
    "bestf(t) = model(a=abest,b=bbest) #Define the best fit function of the form in \"model\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the function bestf(t) with these best choices for a and b.\n",
    "plt2 = plot(bestf(t), t, 1.1, 5.5);\n",
    "pp = plt1 + plt2;\n",
    "show(plt1+plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}