diff --git a/HW6/.~lock.Primary_Energy_monthly.csv# b/HW6/.~lock.Primary_Energy_monthly.csv#
new file mode 100644
index 0000000..eef216f
--- /dev/null
+++ b/HW6/.~lock.Primary_Energy_monthly.csv#
@@ -0,0 +1 @@
+,ryan,fermi,31.03.2017 16:47,file:///home/ryan/.config/libreoffice/4;
\ No newline at end of file
diff --git a/HW6/README.md b/HW6/README.md
new file mode 100644
index 0000000..1f3c86c
--- /dev/null
+++ b/HW6/README.md
@@ -0,0 +1,86 @@
+# Homework #6
+## due 4/14 by 11:59pm
+
+
+0. Create a new github repository called 'curve_fitting'.
+
+ a. Add rcc02007 and pez16103 as collaborators.
+
+ b. Clone the repository to your computer.
+
+
+1. Create a least-squares function called `least_squares.m` that accepts a Z-matrix and
+dependent variable y as input and returns the vector of best-fit constants, a, the
+best-fit function evaluated at each point $f(x_{i})$, and the coefficient of
+determination, r2.
+
+```matlab
+[a,fx,r2]=least_squares(Z,y);
+```
+
+ Test your function on the sets of data in script `problem_1_data.m` and show that the
+ following functions are the best fit lines:
+
+ a. y=0.3745+0.98644x+0.84564/x
+
+ b. y=-11.4887+7.143817x-1.04121 x^2+0.046676 x^3
+
+ c. y=4.0046e^(-1.5x)+2.9213e^(-0.3x)+1.5647e^(-0.05x)
+
+
+2. Use the Temperature and failure data from the Challenger O-rings in lecture_18
+(challenger_oring.csv). Your independent variable is temerature and your dependent
+variable is failure (1=fail, 0=pass). Create a function called `cost_logistic.m` that
+takes the vector `a`, and independent variable `x` and dependent variable `y`. Use the
+function, $\sigma(t)=\frac{1}{1+e^{-t}}$ where $t=a_{0}+a_{1}x$. Use the cost function,
+
+ $J(a_{0},a_{1})=\sum_{i=1}^{n}\left[-y_{i}\log(\sigma(t_{i}))-(1-y_{i})\log((1-\sigma(t_{i})))\right]$
+
+ and gradient
+
+ $\frac{\partial J}{\partial a_{i}}=
+ 1/m\sum_{k=1}^{N}\left(\sigma(t_{k})-y_{k}\right)t_{k}$
+
+ a. edit `cost_logistic.m` so that the output is `[J,grad]` or [cost, gradient]
+
+ b. use the following code to solve for a0 and a1
+
+```matlab
+% Set options for fminunc
+options = optimset('GradObj', 'on', 'MaxIter', 400);
+% Run fminunc to obtain the optimal theta
+% This function will return theta and the cost
+[theta, cost] = ...
+fminunc(@(a)(costFunction(a, x, y)), initial_a, options);
+```
+
+ c. plot the data and the best-fit logistic regression model
+
+```matlab
+plot(x,y, x, sigma(a(1)+a(2)*x))
+```
+
+3. The vertical stress under a corner of a rectangular area subjected to a uniform load of
+intensity $q$ is given by the solution of the Boussinesq's equation:
+
+ $\sigma_{z} =
+ \frac{q}{4\pi}\left(\frac{2mn\sqrt{m^{2}+n^{2}+1}}{m^{2}+n^{2}+1+m^{2}n^{2}}\frac{m^{2}+n^{2}+2}{m^{2}+n^{2}+1}+sin^{-1}\left(\frac{2mn\sqrt{m^{2}+n^{2}+1}}{m^{2}+n^{2}+1+m^{2}n^{2}}\right)\right)$
+
+ Typically, this equation is solved as a table of values where:
+
+ $\sigma_{z}=q f(m,n)$
+
+ where $f(m,n)$ is the influence value, q is the uniform load, m=a/z, n=b/z, a and b are
+ width and length of the rectangular area and z is the depth below the area.
+
+ a. Finish the function `boussinesq_lookup.m` so that when you enter a force, q,
+ dimensions of rectangular area a, b, and depth, z, it uses a third-order polynomial
+ interpolation of the four closest values of m to determine the stress in the vertical
+ direction, sigma_z=$\sigma_{z}$. Use a $0^{th}$-order, polynomial interpolation for
+ the value of n (i.e. round to the closest value of n).
+
+ b. Copy the `boussinesq_lookup.m` to a file called `boussinesq_spline.m` and use a
+ cubic spline to interpolate in two dimensions, both m and n, that returns sigma_z.
+
+
+
diff --git a/HW6/README.pdf b/HW6/README.pdf
new file mode 100644
index 0000000..3cc8991
Binary files /dev/null and b/HW6/README.pdf differ
diff --git a/HW6/boussinesq_lookup.m b/HW6/boussinesq_lookup.m
new file mode 100644
index 0000000..004c91c
--- /dev/null
+++ b/HW6/boussinesq_lookup.m
@@ -0,0 +1,22 @@
+function sigma_z=boussinesq_lookup(q,a,b,z)
+ % function that determines stress under corner of an a by b rectangular platform
+ % z-meters below the platform. The calculated solutions are in the fmn data
+ % m=fmn(:,1)
+ % in column 2, fmn(:,2), n=1.2
+ % in column 3, fmn(:,2), n=1.4
+ % in column 4, fmn(:,2), n=1.6
+
+ fmn= [0.1,0.02926,0.03007,0.03058
+ 0.2,0.05733,0.05894,0.05994
+ 0.3,0.08323,0.08561,0.08709
+ 0.4,0.10631,0.10941,0.11135
+ 0.5,0.12626,0.13003,0.13241
+ 0.6,0.14309,0.14749,0.15027
+ 0.7,0.15703,0.16199,0.16515
+ 0.8,0.16843,0.17389,0.17739];
+
+ m=a/z;
+ n=b/z;
+
+ %...
+end
diff --git a/HW6/cost_logistic.m b/HW6/cost_logistic.m
new file mode 100644
index 0000000..169718f
--- /dev/null
+++ b/HW6/cost_logistic.m
@@ -0,0 +1,27 @@
+function [J, grad] = cost_logistic(a, x, y)
+% cost_logistic Compute cost and gradient for logistic regression
+% J = cost_logistic(theta, X, y) computes the cost of using theta as the
+% parameter for logistic regression and the gradient of the cost
+% w.r.t. to the parameters.
+
+% Initialize some useful values
+N = length(y); % number of training examples
+
+% You need to return the following variables correctly
+J = 0;
+grad = zeros(size(theta));
+
+% ====================== YOUR CODE HERE ======================
+% Instructions: Compute the cost of a particular choice of a.
+% Compute the partial derivatives and set grad to the partial
+% derivatives of the cost w.r.t. each parameter in theta
+%
+% Note: grad should have the same dimensions as theta
+%
+
+
+
+% =============================================================
+
+end
+
diff --git a/HW6/problem_1_data.m b/HW6/problem_1_data.m
new file mode 100644
index 0000000..6af2956
--- /dev/null
+++ b/HW6/problem_1_data.m
@@ -0,0 +1,15 @@
+
+% part a
+xa=[1 2 3 4 5]';
+yb=[2.2 2.8 3.6 4.5 5.5]';
+
+% part b
+
+xb=[3 4 5 7 8 9 11 12]';
+yb=[1.6 3.6 4.4 3.4 2.2 2.8 3.8 4.6]';
+
+% part c
+
+xc=[0.5 1 2 3 4 5 6 7 9];
+yc=[6 4.4 3.2 2.7 2.2 1.9 1.7 1.4 1.1];
+
diff --git a/lecture_16/lecture_16.ipynb b/lecture_16/lecture_16.ipynb
index 7ab08d2..206e84c 100644
--- a/lecture_16/lecture_16.ipynb
+++ b/lecture_16/lecture_16.ipynb
@@ -20644,7 +20644,7 @@
},
{
"cell_type": "code",
- "execution_count": 123,
+ "execution_count": 200,
"metadata": {
"collapsed": false
},
@@ -20793,12 +20793,20 @@
"\n",
"\n",
"\n",
- "\tgnuplot_plot_1a\n",
- "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\tdata\n",
+ "\n",
"\n",
"\n",
- "\t \n",
- "\t\n",
+ "\t\n",
+ "\t\tdata\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\n",
"\t\n",
- "\tgnuplot_plot_2a\n",
+ "\tbest-fit\n",
+ "\n",
+ "\t\n",
+ "\t\tbest-fit\n",
+ "\t\n",
+ "\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
"\n",
"\n",
@@ -20851,6 +20865,7 @@
"a=A\\b\n",
"\n",
"plot(x,y,'o',x,a(1)+a(2)*x)\n",
+ "legend('data','best-fit','Location','NorthWest')\n",
"xlabel('Force (N)')\n",
"ylabel('velocity (m/s)')"
]
diff --git a/lecture_17/.ipynb_checkpoints/lecture_16-checkpoint.ipynb b/lecture_17/.ipynb_checkpoints/lecture_16-checkpoint.ipynb
new file mode 100644
index 0000000..37443ae
--- /dev/null
+++ b/lecture_17/.ipynb_checkpoints/lecture_16-checkpoint.ipynb
@@ -0,0 +1,754 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 171,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 172,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "### Project Ideas so far\n",
+ "\n",
+ "- Nothing yet...probably something heat transfer related\n",
+ "\n",
+ "- Modeling Propulsion or Propagation of Sound Waves\n",
+ "\n",
+ "- Low Thrust Orbital Transfer\n",
+ "\n",
+ "- Tracking motion of a satellite entering orbit until impact\n",
+ "\n",
+ "- What ever you think is best.\n",
+ "\n",
+ "- You had heat transfer project as an option; that sounded cool\n",
+ "\n",
+ "- Heat transfer through a pipe\n",
+ "\n",
+ "- I would prefer to do something with beam/plate mechanics or vibrations than a heat transfer or thermo problem\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Questions from you:\n",
+ "\n",
+ "- Is attempting to divide by zero an acceptable project idea?\n",
+ "\n",
+ "- Would it be alright if we worked in a group of 4?\n",
+ "\n",
+ "- What are acceptable project topics?\n",
+ "\n",
+ "- How do the exams look? \n",
+ "\n",
+ "- Is there no pdf for the lecture today?\n",
+ "\n",
+ "- Thank you for making the formatted lectures available!\n",
+ "\n",
+ "- did you do anything cool over spring break?\n",
+ "\n",
+ "- Could we have a group of 4? We don't want to have to choose which one of us is on their own.\n",
+ "\n",
+ "- In HW 5 should there be 4 vectors as an input?\n",
+ "\n",
+ "- Would it be possible for me to join a group of 3? I seem to be the odd man out in two 3 member groups that my friends are in."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Linear Least Squares Regression\n",
+ "\n",
+ "When approximating a set of data as a polynomial, there will always be error introduced (except in 2 cases). \n",
+ "\n",
+ "For a straight line, the actual data points, $y_{i}$ as a function of the independent variable, $x_{i}$, is:\n",
+ "\n",
+ "$y_{i}=a_{0}+a_{1}x_{i}+e_{i}$\n",
+ "\n",
+ "where $a_{0}$ and $a_{1}$ are the intercept and slope of the line and $e_{i}$ is the error between the approximate function and the recorded data point. \n",
+ "\n",
+ "We make the following assumptions in this analysis:\n",
+ "\n",
+ "1. Each x has a fixed value; it is not random and is known without error.\n",
+ "\n",
+ "2. The y values are independent random variables and all have the same variance.\n",
+ "\n",
+ "3. The y values for a given x must be normally distributed.\n",
+ "\n",
+ "The total error is \n",
+ "\n",
+ "$\\sum_{i=1}^{n}e_{i}=\\sum_{i=1}^{n}(y_{i}-a_{0}-a_{1}x_{i})$\n",
+ "\n",
+ "we don't care about the sign though. One approach has been demonstrated to provide a unique solution is minimizing the sum of squares error or\n",
+ "\n",
+ "$S_{r}=\\sum_{i=1}^{n}e_{i}^{2}=\\sum_{i=1}^{n}(y_{i}-a_{0}-a_{1}x_{i})^{2}$\n",
+ "\n",
+ "Where, $S_{r}$ is the sum of squares error (SSE). \n",
+ "\n",
+ "$\\frac{\\partial S_{r}}{\\partial a_{0}}=-2\\sum(y_{i}-a_{0}-a_{1}x_{i})$\n",
+ "\n",
+ "$\\frac{\\partial S_{r}}{\\partial a_{1}}=-2\\sum(y_{i}-a_{0}-a_{1}x_{i})x_{i}$\n",
+ "\n",
+ "The minimum $S_{r}$ occurrs when the partial derivatives are 0. \n",
+ "\n",
+ "$\\sum y_{i}= \\sum a_{0}+\\sum a_{1}x_{i}$\n",
+ "\n",
+ "$\\sum x_{i}y_{i}= \\sum a_{0}x_{i}+\\sum a_{1}x_{i}^{2}$\n",
+ "\n",
+ "$\\left[\\begin{array}{c}\n",
+ "\\sum y_{i}\\\\\n",
+ "\\sum x_{i}y_{i}\\end{array}\\right]=\n",
+ "\\left[\\begin{array}{cc}\n",
+ "n & \\sum x_{i}\\\\\n",
+ "\\sum x_{i} & \\sum x_{i}^{2}\\end{array}\\right]\n",
+ "\\left[\\begin{array}{c}\n",
+ "a_{0}\\\\\n",
+ "a_{1}\\end{array}\\right]$\n",
+ "\n",
+ "\n",
+ "$b=Ax$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example \n",
+ "\n",
+ "Find drag coefficient with best-fit line to experimental data\n",
+ "\n",
+ "|i | v (m/s) | F (N) |\n",
+ "|---|---|---|\n",
+ "|1 | 10 | 25 |\n",
+ "|2 | 20 | 70 |\n",
+ "|3 | 30 | 380|\n",
+ "|4 | 40 | 550|\n",
+ "|5 | 50 | 610|\n",
+ "|6 | 60 | 1220|\n",
+ "|7 | 70 | 830 |\n",
+ "|8 |80 | 1450|"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 200,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " -234.286\n",
+ " 19.470\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "drag_data=[...\n",
+ "1 , 10 , 25 \n",
+ "2 , 20 , 70 \n",
+ "3 , 30 , 380\n",
+ "4 , 40 , 550\n",
+ "5 , 50 , 610\n",
+ "6 , 60 , 1220\n",
+ "7 , 70 , 830 \n",
+ "8 ,80 , 1450];\n",
+ "x=drag_data(:,2);\n",
+ "y=drag_data(:,3);\n",
+ "\n",
+ "b=[sum(y);sum(x.*y)];\n",
+ "A=[length(x),sum(x);\n",
+ " sum(x), sum(x.^2)];\n",
+ " \n",
+ "a=A\\b\n",
+ "\n",
+ "plot(x,y,'o',x,a(1)+a(2)*x)\n",
+ "legend('data','best-fit','Location','NorthWest')\n",
+ "xlabel('Force (N)')\n",
+ "ylabel('velocity (m/s)')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "How do we know its a \"good\" fit? \n",
+ "\n",
+ "Can compare the sum of squares error to the total sum of squares of the dependent variable (here F). \n",
+ "\n",
+ "$S_{r}=\\sum(y_{i}-a_{0}-a_{1}x_{i})^{2}$\n",
+ "\n",
+ "$S_{t}=\\sum(y_{i}-\\bar{y})^{2}$\n",
+ "\n",
+ "Then, we can calculate the *coefficient of determination*, $r^{2}$ or *correlation coefficient*, r. \n",
+ "\n",
+ "$r^{2}=\\frac{S_{t}-S_{r}}{S_{t}}$\n",
+ "\n",
+ "This represents the relative improvement of assuming that y is a function of x (if the x-values are not random and the y-values are random)\n",
+ "\n",
+ "For further information regarding statistical work on regression, look at \n",
+ "[NIST Statistics Handbook](http://www.itl.nist.gov/div898/handbook/pmd/section4/pmd44.htm)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 128,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "St = 1.8083e+06\n",
+ "St = 1.8083e+06\n"
+ ]
+ }
+ ],
+ "source": [
+ "Sr=sum((y-a(1)-a(2)*x).^2);\n",
+ "St=std(y)^2*(length(y)-1)\n",
+ "St=sum((y-mean(y)).^2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 130,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "r2 = 0.88049\n",
+ "r = 0.93834\n"
+ ]
+ }
+ ],
+ "source": [
+ "r2=(St-Sr)/St\n",
+ "r=sqrt(r2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Limiting cases \n",
+ "\n",
+ "#### $r^{2}=0$ $S_{r}=S{t}$\n",
+ "\n",
+ "#### $r^{2}=1$ $S_{r}=0$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 152,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "gen_examples\n",
+ "plot(x1(1:50:end),y1(1:50:end),'s',x2,y2,'o')\n",
+ "legend('Case 1','Case 2','Location','NorthWest')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 157,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a1 =\n",
+ "\n",
+ " 0.0497269\n",
+ " 0.0016818\n",
+ "\n",
+ "a2 =\n",
+ "\n",
+ " 0\n",
+ " 1\n",
+ "\n",
+ "Sr1 = 0.82607\n",
+ "St1 = 0.82631\n",
+ "coefficient of determination in Case 1 is 0.000286\n",
+ "Sr2 = 0\n",
+ "St2 = 1.0185\n",
+ "coefficient of determination in Case 2 is 1.000000\n"
+ ]
+ }
+ ],
+ "source": [
+ "b1=[sum(y1);sum(x1.*y1)];\n",
+ "A1=[length(x1),sum(x1);\n",
+ " sum(x1), sum(x1.^2)];\n",
+ " \n",
+ "a1=A1\\b1\n",
+ "\n",
+ "b2=[sum(y2);sum(x2.*y2)];\n",
+ "A2=[length(x2),sum(x2);\n",
+ " sum(x2), sum(x2.^2)];\n",
+ " \n",
+ "a2=A2\\b2\n",
+ "\n",
+ "Sr1=sum((y1-a1(1)-a1(2)*x1).^2)\n",
+ "St1=sum((y1-mean(y1)).^2)\n",
+ "fprintf('coefficient of determination in Case 1 is %f\\n',1-Sr1/St1)\n",
+ "\n",
+ "Sr2=sum((y2-a2(1)-a2(2)*x2).^2)\n",
+ "St2=sum((y2-mean(y2)).^2)\n",
+ "\n",
+ "fprintf('coefficient of determination in Case 2 is %f\\n',1-Sr2/St2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "# General Regression (Linear and nonlinear)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb b/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb
new file mode 100644
index 0000000..2fd6442
--- /dev/null
+++ b/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_17/gen_examples.m b/lecture_17/gen_examples.m
new file mode 100644
index 0000000..c5d6890
--- /dev/null
+++ b/lecture_17/gen_examples.m
@@ -0,0 +1,5 @@
+x1=linspace(0,1,1000)';
+y1=0*x1+0.1*rand(length(x1),1);
+
+x2=linspace(0,1,10)';
+y2=1*x2;
diff --git a/lecture_17/in-class_regression.pdf b/lecture_17/in-class_regression.pdf
new file mode 100644
index 0000000..628b495
Binary files /dev/null and b/lecture_17/in-class_regression.pdf differ
diff --git a/lecture_17/lecture_17.ipynb b/lecture_17/lecture_17.ipynb
new file mode 100644
index 0000000..35e61da
--- /dev/null
+++ b/lecture_17/lecture_17.ipynb
@@ -0,0 +1,2710 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "### Project Ideas so far\n",
+ "\n",
+ "- Nothing yet...probably something heat transfer related\n",
+ "\n",
+ "- Modeling Propulsion or Propagation of Sound Waves\n",
+ "\n",
+ "- Low Thrust Orbital Transfer\n",
+ "\n",
+ "- Tracking motion of a satellite entering orbit until impact\n",
+ "\n",
+ "- What ever you think is best.\n",
+ "\n",
+ "- You had heat transfer project as an option; that sounded cool\n",
+ "\n",
+ "- Heat transfer through a pipe\n",
+ "\n",
+ "- I would prefer to do something with beam/plate mechanics or vibrations than a heat transfer or thermo problem\n",
+ "\n",
+ "- Modeling Couette flow with a pressure gradient using a discretized form of the Navier-Stokes equation and comparing it to the analytical solution\n",
+ "\n",
+ "- Software to instruct a robotic arm to orient itself based on input from a gyroscope"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Questions from you:\n",
+ "\n",
+ "How was your spring break\n",
+ "\n",
+ "Are grades for hw 3 and 4 going to be posted?\n",
+ "\n",
+ "For class occasionally switching to doc cam and working through problems by hand might help get ideas across.\n",
+ "\n",
+ "Are you assigning those who do not have groups to groups?\n",
+ "\n",
+ "How do you graph a standard distribution in Matlab?\n",
+ "\n",
+ "What is the longest code you have written?\n",
+ "\n",
+ "how to approach probability for hw 5?\n",
+ "??\n",
+ "\n",
+ "Will you be assigning groups to people who do not currently have one? \n",
+ "\n",
+ "what are some basic guidelines we should use to brainstorm project ideas?\n",
+ "\n",
+ "Are you a fan of bananas?\n",
+ "\n",
+ "Going through code isn't the most helpful, because you can easily lose interest. But I am not sure what else you can do.\n",
+ "\n",
+ "Has lecture 15 been posted yet? Looking in the repository I can't seem to find it."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Linear Least Squares Regression\n",
+ "\n",
+ "When approximating a set of data as a polynomial, there will always be error introduced (except in 2 cases). \n",
+ "\n",
+ "For a straight line, the actual data points, $y_{i}$ as a function of the independent variable, $x_{i}$, is:\n",
+ "\n",
+ "$y_{i}=a_{0}+a_{1}x_{i}+e_{i}$\n",
+ "\n",
+ "where $a_{0}$ and $a_{1}$ are the intercept and slope of the line and $e_{i}$ is the error between the approximate function and the recorded data point. \n",
+ "\n",
+ "We make the following assumptions in this analysis:\n",
+ "\n",
+ "1. Each x has a fixed value; it is not random and is known without error.\n",
+ "\n",
+ "2. The y values are independent random variables and all have the same variance.\n",
+ "\n",
+ "3. The y values for a given x must be normally distributed.\n",
+ "\n",
+ "The total error is \n",
+ "\n",
+ "$\\sum_{i=1}^{n}e_{i}=\\sum_{i=1}^{n}(y_{i}-a_{0}-a_{1}x_{i})$\n",
+ "\n",
+ "we don't care about the sign though. One approach has been demonstrated to provide a unique solution is minimizing the sum of squares error or\n",
+ "\n",
+ "$S_{r}=\\sum_{i=1}^{n}e_{i}^{2}=\\sum_{i=1}^{n}(y_{i}-a_{0}-a_{1}x_{i})^{2}$\n",
+ "\n",
+ "Where, $S_{r}$ is the sum of squares error (SSE). \n",
+ "\n",
+ "$\\frac{\\partial S_{r}}{\\partial a_{0}}=-2\\sum(y_{i}-a_{0}-a_{1}x_{i})$\n",
+ "\n",
+ "$\\frac{\\partial S_{r}}{\\partial a_{1}}=-2\\sum(y_{i}-a_{0}-a_{1}x_{i})x_{i}$\n",
+ "\n",
+ "The minimum $S_{r}$ occurrs when the partial derivatives are 0. \n",
+ "\n",
+ "$\\sum y_{i}= \\sum a_{0}+\\sum a_{1}x_{i}$\n",
+ "\n",
+ "$\\sum x_{i}y_{i}= \\sum a_{0}x_{i}+\\sum a_{1}x_{i}^{2}$\n",
+ "\n",
+ "$\\left[\\begin{array}{c}\n",
+ "\\sum y_{i}\\\\\n",
+ "\\sum x_{i}y_{i}\\end{array}\\right]=\n",
+ "\\left[\\begin{array}{cc}\n",
+ "n & \\sum x_{i}\\\\\n",
+ "\\sum x_{i} & \\sum x_{i}^{2}\\end{array}\\right]\n",
+ "\\left[\\begin{array}{c}\n",
+ "a_{0}\\\\\n",
+ "a_{1}\\end{array}\\right]$\n",
+ "\n",
+ "\n",
+ "$b=Ax$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example \n",
+ "\n",
+ "Find drag coefficient with best-fit line to experimental data\n",
+ "\n",
+ "|i | v (m/s) | F (N) |\n",
+ "|---|---|---|\n",
+ "|1 | 10 | 25 |\n",
+ "|2 | 20 | 70 |\n",
+ "|3 | 30 | 380|\n",
+ "|4 | 40 | 550|\n",
+ "|5 | 50 | 610|\n",
+ "|6 | 60 | 1220|\n",
+ "|7 | 70 | 830 |\n",
+ "|8 |80 | 1450|"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 68,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " -234.286\n",
+ " 19.470\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "drag_data=[...\n",
+ "1 , 10 , 25 \n",
+ "2 , 20 , 70 \n",
+ "3 , 30 , 380\n",
+ "4 , 40 , 550\n",
+ "5 , 50 , 610\n",
+ "6 , 60 , 1220\n",
+ "7 , 70 , 830 \n",
+ "8 ,80 , 1450];\n",
+ "x=drag_data(:,2);\n",
+ "y=drag_data(:,3);\n",
+ "\n",
+ "b=[sum(y);sum(x.*y)];\n",
+ "A=[length(x),sum(x);\n",
+ " sum(x), sum(x.^2)];\n",
+ " \n",
+ "a=A\\b\n",
+ "\n",
+ "plot(x,y,'o',x,a(1)+a(2)*x)\n",
+ "legend('data','best-fit','Location','NorthWest')\n",
+ "xlabel('Force (N)')\n",
+ "ylabel('velocity (m/s)')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 72,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(x,y-a(1)-40*x)\n",
+ "legend('data','best-fit','Location','NorthWest')\n",
+ "ylabel('residuals (N)')\n",
+ "xlabel('velocity (m/s)')\n",
+ "title('Model does not capture measurements')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "How do we know its a \"good\" fit? \n",
+ "\n",
+ "Can compare the sum of squares error to the total sum of squares of the dependent variable (here F). \n",
+ "\n",
+ "$S_{r}=\\sum(y_{i}-a_{0}-a_{1}x_{i})^{2}$\n",
+ "\n",
+ "$S_{t}=\\sum(y_{i}-\\bar{y})^{2}$\n",
+ "\n",
+ "Then, we can calculate the *coefficient of determination*, $r^{2}$ or *correlation coefficient*, r. \n",
+ "\n",
+ "$r^{2}=\\frac{S_{t}-S_{r}}{S_{t}}$\n",
+ "\n",
+ "This represents the relative improvement of assuming that y is a function of x (if the x-values are not random and the y-values are random)\n",
+ "\n",
+ "For further information regarding statistical work on regression, look at \n",
+ "[NIST Statistics Handbook](http://www.itl.nist.gov/div898/handbook/pmd/section4/pmd44.htm)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 73,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "St = 1.8083e+06\n",
+ "St = 1.8083e+06\n"
+ ]
+ }
+ ],
+ "source": [
+ "Sr=sum((y-a(1)-a(2)*x).^2);\n",
+ "St=std(y)^2*(length(y)-1)\n",
+ "St=sum((y-mean(y)).^2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 74,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "r2 = 0.88049\n",
+ "r = 0.93834\n"
+ ]
+ }
+ ],
+ "source": [
+ "r2=(St-Sr)/St\n",
+ "r=sqrt(r2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Limiting cases \n",
+ "\n",
+ "#### $r^{2}=0$ $S_{r}=S{t}$\n",
+ "\n",
+ "#### $r^{2}=1$ $S_{r}=0$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 76,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "gen_examples\n",
+ "plot(x1(1:50:end),y1(1:50:end),'s',x2,y2,'o')\n",
+ "legend('Case 1','Case 2','Location','NorthWest')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 77,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a1 =\n",
+ "\n",
+ " 0.0482496\n",
+ " 0.0017062\n",
+ "\n",
+ "a2 =\n",
+ "\n",
+ " 0\n",
+ " 1\n",
+ "\n",
+ "Sr1 = 0.83505\n",
+ "St1 = 0.83529\n",
+ "coefficient of determination in Case 1 is 0.000291\n",
+ "Sr2 = 0\n",
+ "St2 = 1.0185\n",
+ "coefficient of determination in Case 2 is 1.000000\n"
+ ]
+ }
+ ],
+ "source": [
+ "b1=[sum(y1);sum(x1.*y1)];\n",
+ "A1=[length(x1),sum(x1);\n",
+ " sum(x1), sum(x1.^2)];\n",
+ " \n",
+ "a1=A1\\b1\n",
+ "\n",
+ "b2=[sum(y2);sum(x2.*y2)];\n",
+ "A2=[length(x2),sum(x2);\n",
+ " sum(x2), sum(x2.^2)];\n",
+ " \n",
+ "a2=A2\\b2\n",
+ "\n",
+ "Sr1=sum((y1-a1(1)-a1(2)*x1).^2)\n",
+ "St1=sum((y1-mean(y1)).^2)\n",
+ "fprintf('coefficient of determination in Case 1 is %f\\n',1-Sr1/St1)\n",
+ "\n",
+ "Sr2=sum((y2-a2(1)-a2(2)*x2).^2)\n",
+ "St2=sum((y2-mean(y2)).^2)\n",
+ "\n",
+ "fprintf('coefficient of determination in Case 2 is %f\\n',1-Sr2/St2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "# General Linear Regression"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In general, you may want to fit other polynomials besides degree-1 (straight-lines)\n",
+ "\n",
+ "$y=a_{0}+a_{1}x+a_{2}x^{2}+\\cdots+a_{m}x^{m}+e$\n",
+ "\n",
+ "Now, the solution for $a_{0},~a_{1},...a_{m}$ is the minimization of m+1-dependent linear equations. \n",
+ "\n",
+ "Consider the following data:\n",
+ "\n",
+ "| x | y |\n",
+ "|---|---|\n",
+ "| 0.00 | 21.50 |\n",
+ "| 2.00 | 20.84 |\n",
+ "| 4.00 | 23.19 |\n",
+ "| 6.00 | 22.69 |\n",
+ "| 8.00 | 30.27 |\n",
+ "| 10.00 | 40.11 |\n",
+ "| 12.00 | 43.31 |\n",
+ "| 14.00 | 54.79 |\n",
+ "| 16.00 | 70.88 |\n",
+ "| 18.00 | 89.48 |\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 78,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "load xy_data.csv\n",
+ "x=xy_data(1,:)';\n",
+ "y=xy_data(2,:)';\n",
+ "plot(x,y,'o')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 46,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "xy_data =\n",
+ "\n",
+ " Columns 1 through 7:\n",
+ "\n",
+ " 0.00000 2.00000 4.00000 6.00000 8.00000 10.00000 12.00000\n",
+ " 21.50114 20.84153 23.19201 22.69498 30.26687 40.11075 43.30543\n",
+ "\n",
+ " Columns 8 through 11:\n",
+ "\n",
+ " 14.00000 16.00000 18.00000 20.00000\n",
+ " 54.78730 70.88443 89.48368 97.28135\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "xy_data"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The model can be rewritten as \n",
+ "\n",
+ "$y=\\left[Z\\right]a+e$\n",
+ "\n",
+ "where $a=\\left[\\begin{array}{c}\n",
+ " a_{0}\\\\\n",
+ " a_{1}\\\\\n",
+ " a_{2}\\end{array}\\right]$\n",
+ " \n",
+ "$[Z]=\\left[\\begin{array} \n",
+ "1 & x_{1} & x_{1}^{2} \\\\\n",
+ "1 & x_{2} & x_{2}^{2} \\\\\n",
+ "1 & x_{3} & x_{3}^{2} \\\\\n",
+ "1 & x_{4} & x_{4}^{2} \\\\\n",
+ "1 & x_{5} & x_{5}^{2} \\\\\n",
+ "1 & x_{6} & x_{6}^{2} \\\\\n",
+ "1 & x_{7} & x_{7}^{2} \\\\\n",
+ "1 & x_{8} & x_{8}^{2} \\\\\n",
+ "1 & x_{9} & x_{9}^{2} \\\\\n",
+ "1 & x_{10} & x_{10}^{2} \\end{array}\\right]$\n",
+ "\n",
+ "The sum of squares residuals for this model is\n",
+ "\n",
+ "$S_{r}=\\sum_{i=1}^{n}\\left(y_{i}-\\sum_{j=0}^{m}a_{j}z_{ji}\\right)$\n",
+ "\n",
+ "Minimizing this function results in\n",
+ "\n",
+ "$y=[Z]a$\n",
+ "\n",
+ "->**A standard Linear Algebra Problem**\n",
+ "\n",
+ "*the vector a is unknown, and Z is calculated based upon the assumed function*\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 79,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Z =\n",
+ "\n",
+ " 1 0 0\n",
+ " 1 2 4\n",
+ " 1 4 16\n",
+ " 1 6 36\n",
+ " 1 8 64\n",
+ " 1 10 100\n",
+ " 1 12 144\n",
+ " 1 14 196\n",
+ " 1 16 256\n",
+ " 1 18 324\n",
+ " 1 20 400\n",
+ "\n",
+ "a =\n",
+ "\n",
+ " 21.40341\n",
+ " -0.81538\n",
+ " 0.23935\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "Z=[x.^0,x,x.^2]\n",
+ "\n",
+ "a=Z\\y"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 80,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x_fcn=linspace(min(x),max(x));\n",
+ "plot(x,y,'o',x_fcn,a(1)+a(2)*x_fcn+a(3)*x_fcn.^2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### General Coefficient of Determination\n",
+ "\n",
+ "$r^{2}=\\frac{S_{t}-S_{r}}{S_{t}}=1-\\frac{S_{r}}{S_t}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 49,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "St = 27.923\n",
+ "Sr = 2.6366\n"
+ ]
+ }
+ ],
+ "source": [
+ "St=std(y)\n",
+ "Sr=std(y-a(1)-a(2)*x-a(3)*x.^2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 81,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "the coefficient of determination for this fit is 0.880485\n",
+ "the correlation coefficient this fit is 0.938342\n"
+ ]
+ }
+ ],
+ "source": [
+ "r2=1-Sr/St;\n",
+ "r=sqrt(r2);\n",
+ "\n",
+ "fprintf('the coefficient of determination for this fit is %f',r2)\n",
+ "fprintf('the correlation coefficient this fit is %f',r)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Compare this to a straight line fit"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 83,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "St = 27.923\n",
+ "Sr = 9.2520\n",
+ "the coefficient of determination for this fit is 0.668655\n",
+ "the correlation coefficient this fit is 0.817713\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Z=[ones(size(x)) x];\n",
+ "a_line=Z\\y;\n",
+ "x_fcn=linspace(min(x),max(x));\n",
+ "plot(x,y,'o',x_fcn,a(1)+a(2)*x_fcn+a(3)*x_fcn.^2,...\n",
+ "x_fcn,a_line(1)+a_line(2)*x_fcn)\n",
+ "St=std(y)\n",
+ "Sr=std(y-a_line(1)-a_line(2)*x)\n",
+ "r2=1-Sr/St;\n",
+ "r=sqrt(r2);\n",
+ "\n",
+ "fprintf('the coefficient of determination for this fit is %f',r2)\n",
+ "fprintf('the correlation coefficient this fit is %f',r)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 85,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(x,y-a_line(1)-a_line(2)*x)\n",
+ "title('residuals of straight line')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 86,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(x,y-a(1)-a(2)*x-a(3)*x.^2)\n",
+ "title('residuals of parabola')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Warning \n",
+ "**Coefficient of determination reduction does not always mean a better fit**\n",
+ "\n",
+ "Try the function, \n",
+ "\n",
+ "$y(x)=a_0+a_{1}x+a_{2}x^{2}+a_{4}x^{4}+a_{5}x^{5}+a_{5}x^{5}+a_{6}x^{6}+a_{7}x^{7}+a_{8}x^{8}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 90,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a_overfit =\n",
+ "\n",
+ " 2.1487e+01\n",
+ " -1.4264e+01\n",
+ " 1.5240e+01\n",
+ " -6.0483e+00\n",
+ " 1.1887e+00\n",
+ " -1.2651e-01\n",
+ " 7.4379e-03\n",
+ " -2.2702e-04\n",
+ " 2.8063e-06\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "Z=[ones(size(x)) x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7 x.^8];\n",
+ "a_overfit=Z\\y"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 91,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(x,y,'o',x_fcn,a(1)+a(2)*x_fcn+a(3)*x_fcn.^2,...\n",
+ "x_fcn,a_overfit(1)+a_overfit(2)*x_fcn+a_overfit(3)*x_fcn.^2+...\n",
+ "a_overfit(4)*x_fcn.^3+a_overfit(5)*x_fcn.^4+...\n",
+ "a_overfit(6)*x_fcn.^5+a_overfit(7)*x_fcn.^6+...\n",
+ "a_overfit(8)*x_fcn.^7+a_overfit(9)*x_fcn.^8)\n",
+ "legend('data','parabola','n=8-fit','Location','NorthWest')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 92,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "St = 27.923\n",
+ "Sr = 0.77320\n",
+ "the coefficient of determination for this fit is 0.972309\n",
+ "the correlation coefficient this fit is 0.986057\n"
+ ]
+ }
+ ],
+ "source": [
+ "St=std(y)\n",
+ "Sr=std(y-polyval(a_overfit(end:-1:1),x))\n",
+ "r2=1-Sr/St;\n",
+ "r=sqrt(r2);\n",
+ "\n",
+ "fprintf('the coefficient of determination for this fit is %f',r2)\n",
+ "fprintf('the correlation coefficient this fit is %f',r)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Linear Regression is only limited by the ability to separate the parameters from the function to achieve\n",
+ "\n",
+ "$y=[Z]a$\n",
+ "\n",
+ "$Z$ can be any function of the independent variable(s)\n",
+ "\n",
+ "**Example**:\n",
+ "We assume two functions are added together, sin(t) and sin(3t). What are the amplitudes?\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "load sin_data.csv\n",
+ "t=sin_data(1,:)';\n",
+ "y=sin_data(2,:)';\n",
+ "plot(t,y)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " 0.99727\n",
+ " 0.50251\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "Z=[sin(t) sin(3*t)];\n",
+ "a=Z\\y"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(t,y,'.',t,a(1)*sin(t)+a(2)*sin(3*t))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_17/octave-workspace b/lecture_17/octave-workspace
new file mode 100644
index 0000000..c9ec2aa
Binary files /dev/null and b/lecture_17/octave-workspace differ
diff --git a/lecture_17/q1.png b/lecture_17/q1.png
new file mode 100644
index 0000000..7c89cbf
Binary files /dev/null and b/lecture_17/q1.png differ
diff --git a/lecture_17/q2.png b/lecture_17/q2.png
new file mode 100644
index 0000000..7eea25c
Binary files /dev/null and b/lecture_17/q2.png differ
diff --git a/lecture_17/sin_data.csv b/lecture_17/sin_data.csv
new file mode 100644
index 0000000..a5ff272
--- /dev/null
+++ b/lecture_17/sin_data.csv
@@ -0,0 +1,2 @@
+ 0.00000000e+00 1.26933037e-01 2.53866073e-01 3.80799110e-01 5.07732146e-01 6.34665183e-01 7.61598219e-01 8.88531256e-01 1.01546429e+00 1.14239733e+00 1.26933037e+00 1.39626340e+00 1.52319644e+00 1.65012947e+00 1.77706251e+00 1.90399555e+00 2.03092858e+00 2.15786162e+00 2.28479466e+00 2.41172769e+00 2.53866073e+00 2.66559377e+00 2.79252680e+00 2.91945984e+00 3.04639288e+00 3.17332591e+00 3.30025895e+00 3.42719199e+00 3.55412502e+00 3.68105806e+00 3.80799110e+00 3.93492413e+00 4.06185717e+00 4.18879020e+00 4.31572324e+00 4.44265628e+00 4.56958931e+00 4.69652235e+00 4.82345539e+00 4.95038842e+00 5.07732146e+00 5.20425450e+00 5.33118753e+00 5.45812057e+00 5.58505361e+00 5.71198664e+00 5.83891968e+00 5.96585272e+00 6.09278575e+00 6.21971879e+00 6.34665183e+00 6.47358486e+00 6.60051790e+00 6.72745093e+00 6.85438397e+00 6.98131701e+00 7.10825004e+00 7.23518308e+00 7.36211612e+00 7.48904915e+00 7.61598219e+00 7.74291523e+00 7.86984826e+00 7.99678130e+00 8.12371434e+00 8.25064737e+00 8.37758041e+00 8.50451345e+00 8.63144648e+00 8.75837952e+00 8.88531256e+00 9.01224559e+00 9.13917863e+00 9.26611167e+00 9.39304470e+00 9.51997774e+00 9.64691077e+00 9.77384381e+00 9.90077685e+00 1.00277099e+01 1.01546429e+01 1.02815760e+01 1.04085090e+01 1.05354420e+01 1.06623751e+01 1.07893081e+01 1.09162411e+01 1.10431742e+01 1.11701072e+01 1.12970402e+01 1.14239733e+01 1.15509063e+01 1.16778394e+01 1.18047724e+01 1.19317054e+01 1.20586385e+01 1.21855715e+01 1.23125045e+01 1.24394376e+01 1.25663706e+01
+ 9.15756288e-02 3.39393873e-01 6.28875306e-01 7.67713096e-01 1.05094584e+00 9.70887288e-01 9.84265740e-01 1.02589034e+00 8.53218113e-01 6.90197665e-01 5.51277193e-01 5.01564914e-01 5.25455797e-01 5.87052838e-01 5.41394658e-01 7.12365594e-01 8.14839678e-01 9.80181855e-01 9.44430709e-01 1.06728057e+00 1.15166322e+00 8.99464065e-01 7.77225453e-01 5.92618124e-01 3.08822183e-01 -1.07884730e-03 -3.46563271e-01 -5.64836023e-01 -8.11931510e-01 -1.05925186e+00 -1.13323611e+00 -1.11986890e+00 -8.88336727e-01 -9.54113139e-01 -6.81378679e-01 -6.02369117e-01 -4.78684439e-01 -5.88160325e-01 -4.93580777e-01 -5.68747320e-01 -7.51641934e-01 -8.14672884e-01 -9.53191554e-01 -9.55337518e-01 -9.85995556e-01 -9.63373597e-01 -1.01511061e+00 -7.56467517e-01 -4.17379564e-01 -1.22340361e-01 2.16273929e-01 5.16909714e-01 7.77031694e-01 1.00653798e+00 9.35718089e-01 1.00660116e+00 1.11177057e+00 9.85485116e-01 8.54344900e-01 6.26444042e-01 6.28124048e-01 4.27764254e-01 5.93991751e-01 4.79248018e-01 7.17522492e-01 7.35927848e-01 9.08802925e-01 9.38646871e-01 1.13125860e+00 1.07247935e+00 1.05198782e+00 9.41647332e-01 6.98801244e-01 4.03193543e-01 1.37009682e-01 -1.43203880e-01 -4.64369445e-01 -6.94978252e-01 -1.03483196e+00 -1.10261288e+00 -1.12892727e+00 -1.03902484e+00 -8.53573083e-01 -7.01815315e-01 -6.84745997e-01 -6.14189417e-01 -4.70090797e-01 -5.95052432e-01 -5.96497000e-01 -5.66861911e-01 -7.18239679e-01 -9.52873043e-01 -9.37512847e-01 -1.15782985e+00 -1.03858206e+00 -1.03182712e+00 -8.45121554e-01 -5.61821980e-01 -2.83427014e-01 -8.27056140e-02
diff --git a/lecture_17/xy_data.csv b/lecture_17/xy_data.csv
new file mode 100644
index 0000000..15b919d
--- /dev/null
+++ b/lecture_17/xy_data.csv
@@ -0,0 +1,2 @@
+ 0.00000000e+00 2.00000000e+00 4.00000000e+00 6.00000000e+00 8.00000000e+00 1.00000000e+01 1.20000000e+01 1.40000000e+01 1.60000000e+01 1.80000000e+01 2.00000000e+01
+ 2.15011412e+01 2.08415256e+01 2.31920098e+01 2.26949773e+01 3.02668745e+01 4.01107461e+01 4.33054255e+01 5.47873036e+01 7.08844291e+01 8.94836828e+01 9.72813533e+01
diff --git a/lecture_18/.Newtint.m.swp b/lecture_18/.Newtint.m.swp
new file mode 100644
index 0000000..769fe2b
Binary files /dev/null and b/lecture_18/.Newtint.m.swp differ
diff --git a/lecture_18/.ipynb_checkpoints/lecture_18-checkpoint.ipynb b/lecture_18/.ipynb_checkpoints/lecture_18-checkpoint.ipynb
new file mode 100644
index 0000000..2fd6442
--- /dev/null
+++ b/lecture_18/.ipynb_checkpoints/lecture_18-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_18/Newtint.m b/lecture_18/Newtint.m
new file mode 100644
index 0000000..e4c6c83
--- /dev/null
+++ b/lecture_18/Newtint.m
@@ -0,0 +1,34 @@
+function yint = Newtint(x,y,xx)
+% Newtint: Newton interpolating polynomial
+% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton
+% interpolating polynomial based on n data points (x, y)
+% to determine a value of the dependent variable (yint)
+% at a given value of the independent variable, xx.
+% input:
+% x = independent variable
+% y = dependent variable
+% xx = value of independent variable at which
+% interpolation is calculated
+% output:
+% yint = interpolated value of dependent variable
+
+% compute the finite divided differences in the form of a
+% difference table
+n = length(x);
+if length(y)~=n, error('x and y must be same length'); end
+b = zeros(n,n);
+% assign dependent variables to the first column of b.
+b(:,1) = y(:); % the (:) ensures that y is a column vector.
+for j = 2:n
+ for i = 1:n-j+1
+ b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
+ end
+end
+%b
+% use the finite divided differences to interpolate
+xt = 1;
+yint = b(1,1);
+for j = 1:n-1
+ xt = xt*(xx-x(j));
+ yint = yint+b(1,j+1)*xt;
+end
diff --git a/lecture_18/challenger_oring.csv b/lecture_18/challenger_oring.csv
new file mode 100644
index 0000000..11d647e
--- /dev/null
+++ b/lecture_18/challenger_oring.csv
@@ -0,0 +1,24 @@
+Flight#,Temp,O-Ring Problem
+1,53,1
+2,57,1
+3,58,1
+4,63,1
+5,66,0
+6,66.8,0
+7,67,0
+8,67.2,0
+9,68,0
+10,69,0
+11,69.8,1
+12,69.8,0
+13,70.2,1
+14,70.2,0
+15,72,0
+16,73,0
+17,75,0
+18,75,1
+19,75.8,0
+20,76.2,0
+21,78,0
+22,79,0
+23,81,0
diff --git a/lecture_18/lecture_18.ipynb b/lecture_18/lecture_18.ipynb
new file mode 100644
index 0000000..f1fd1b3
--- /dev/null
+++ b/lecture_18/lecture_18.ipynb
@@ -0,0 +1,2197 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 115,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.081014\n",
+ " 0.317493\n",
+ " 0.690279\n",
+ " 1.169170\n",
+ " 1.715370\n",
+ " 2.284630\n",
+ " 2.830830\n",
+ " 3.309721\n",
+ " 3.682507\n",
+ " 3.918986\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "N=10;\n",
+ "A_beam=diag(ones(N,1))*2+diag(ones(N-1,1)*-1,-1)+diag(ones(N-1,1)*-1,1);\n",
+ "[v,e]=eig(A_beam);\n",
+ "diag(e)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Nonlinear Regression\n",
+ "\n",
+ "We can define any function and minimize the sum of squares error even if the constants cannot be separated.\n",
+ "\n",
+ "$S_{r}=\\left[y-f(z_{1},z_{2},...)\\right]^{2}$\n",
+ "\n",
+ "Consider the function, \n",
+ "\n",
+ "$f(x) = a_{0}(1-e^{a_{1}x})$\n",
+ "\n",
+ "We can define the sum of squares error as a function of $a_{0}$ and $a_{1}$:\n",
+ "\n",
+ "$f_{SSE}(a_{0},a_{1})=\\sum_{i=1}^{n}\\left[y- a_{0}(1-e^{a_{1}x})\\right]^{2}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "function [SSE,yhat] = sse_nonlin_exp(a,x,y)\n",
+ " % This is a sum of squares error function based on \n",
+ " % the two input constants a0 and a1 where a=[a0,a1]\n",
+ " % and the data is x (independent), y (dependent)\n",
+ " % and yhat is the model with the given a0 and a1 values\n",
+ " a0=a(1);\n",
+ " a1=a(2);\n",
+ " yhat=a0*(1-exp(a1*x));\n",
+ " SSE=sum((y-a0*(1-exp(a1*x))).^2);\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Where the data we are fitting is:\n",
+ "\n",
+ "| x | y |\n",
+ "|---|---|\n",
+ " | 0.0 | 0.41213|\n",
+ " | 1.0 | -2.65190|\n",
+ " | 2.0 | 15.04049|\n",
+ " | 3.0 | 5.19368|\n",
+ " | 4.0 | -0.71086|\n",
+ " | 5.0 | 12.69008|\n",
+ " | 6.0 | 29.20309|\n",
+ " | 7.0 | 58.68879|\n",
+ " | 8.0 | 91.61117|\n",
+ " | 9.0 | 173.75492|\n",
+ " | 10.0 | 259.04083|"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "data=[\n",
+ " 0.00000 0.41213\n",
+ " 1.00000 -2.65190\n",
+ " 2.00000 15.04049\n",
+ " 3.00000 5.19368\n",
+ " 4.00000 -0.71086\n",
+ " 5.00000 12.69008\n",
+ " 6.00000 29.20309\n",
+ " 7.00000 58.68879\n",
+ " 8.00000 91.61117\n",
+ " 9.00000 173.75492\n",
+ " 10.00000 259.04083\n",
+ "\n",
+ "];\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 116,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "the sum of squares for a0=-2.00 and a1=0.20 is 98118.4\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[SSE,yhat]=sse_nonlin_exp([-2,0.2],data(:,1),data(:,2));\n",
+ "fprintf('the sum of squares for a0=%1.2f and a1=%1.2f is %1.1f',...\n",
+ "-2,0.2,SSE)\n",
+ "plot(data(:,1),data(:,2),'o',data(:,1),yhat)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 127,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " -1.71891 0.50449\n",
+ "\n",
+ "fsse = 633.70\n"
+ ]
+ }
+ ],
+ "source": [
+ "[a,fsse]=fminsearch(@(a) sse_nonlin_exp(a,data(:,1),data(:,2)),[-2,0.2])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 128,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[sse,yhat]=sse_nonlin_exp(a,data(:,1),data(:,2));\n",
+ "plot(data(:,1),data(:,2),'o',data(:,1),yhat)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 130,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[sse,yhat]=sse_nonlin_exp(a,data(:,1),data(:,2));\n",
+ "plot(data(:,1),data(:,2)-yhat)\n",
+ "title('residuals of function')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Case Study: Logistic Regression\n",
+ "\n",
+ "Many times the variable you predict is a binary (or discrete) value, such as pass/fail, broken/not-broken, etc. \n",
+ "\n",
+ "One method to fit this type of data is called [**logistic regression**](https://en.wikipedia.org/wiki/Logistic_regression).\n",
+ "\n",
+ "[Logistic Regression link 2](http://www.holehouse.org/mlclass/06_Logistic_Regression.html)\n",
+ "\n",
+ "We use a function that varies from 0 to 1 called a logistic function:\n",
+ "\n",
+ "$\\sigma(t)=\\frac{1}{1+e^{-t}}$\n",
+ "\n",
+ "We can use this function to describe the likelihood of failure (1) or success (0). When t=0, the probability of failure is 50%. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 131,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "t=linspace(-10,10);\n",
+ "sigma=@(t) 1./(1+exp(-t));\n",
+ "plot(t,sigma(t))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we make the assumption that we can predict the boundary between the pass-fail criteria with a function of our independent variable e.g.\n",
+ "\n",
+ "$y=\\left\\{\\begin{array}{cc} \n",
+ "1 & a_{0}+a_{1}x +\\epsilon >0 \\\\\n",
+ "0 & else \\end{array} \\right\\}$\n",
+ "\n",
+ "so the logistic function is now:\n",
+ "\n",
+ "$\\sigma(x)=\\frac{1}{1+e^{-(a_{0}+a_{1}x)}}$\n",
+ "\n",
+ "Here, there is not a direct sum of squares error, so we minimize a cost function: \n",
+ "\n",
+ "$J(a_{0},a_{1})=\\sum_{i=1}^{n}\\left[-y_{i}\\log(\\sigma(x_{i}))-(1-y_{i})\\log((1-\\sigma(x_{i})))\\right]$\n",
+ "\n",
+ "y=0,1 \n",
+ "\n",
+ "So the cost function either sums the "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example: Challenger O-ring failures\n",
+ "\n",
+ "The O-rings on the Challenger shuttles had problems when temperatures became low. We can look at the conditions when damage was observed to determine the likelihood of failure. \n",
+ "\n",
+ "[Challenger O-ring data powerpoint](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjZvL7jkP3SAhUp04MKHXkXDkMQFggcMAA&url=http%3A%2F%2Fwww.stat.ufl.edu%2F~winner%2Fcases%2Fchallenger.ppt&usg=AFQjCNFyjwT7NmRthDkDEgch75Fc5dc66w&sig2=_qeteX6-ZEBwPW8SZN1mIA)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 132,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "oring =\n",
+ "\n",
+ " 1.00000 53.00000 1.00000\n",
+ " 2.00000 57.00000 1.00000\n",
+ " 3.00000 58.00000 1.00000\n",
+ " 4.00000 63.00000 1.00000\n",
+ " 5.00000 66.00000 0.00000\n",
+ " 6.00000 66.80000 0.00000\n",
+ " 7.00000 67.00000 0.00000\n",
+ " 8.00000 67.20000 0.00000\n",
+ " 9.00000 68.00000 0.00000\n",
+ " 10.00000 69.00000 0.00000\n",
+ " 11.00000 69.80000 1.00000\n",
+ " 12.00000 69.80000 0.00000\n",
+ " 13.00000 70.20000 1.00000\n",
+ " 14.00000 70.20000 0.00000\n",
+ " 15.00000 72.00000 0.00000\n",
+ " 16.00000 73.00000 0.00000\n",
+ " 17.00000 75.00000 0.00000\n",
+ " 18.00000 75.00000 1.00000\n",
+ " 19.00000 75.80000 0.00000\n",
+ " 20.00000 76.20000 0.00000\n",
+ " 21.00000 78.00000 0.00000\n",
+ " 22.00000 79.00000 0.00000\n",
+ " 23.00000 81.00000 0.00000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "% read data from csv file \n",
+ "% col 1 = index\n",
+ "% col 2 = temperature\n",
+ "% col 3 = 1 if damaged, 0 if undamaged\n",
+ "oring=dlmread('challenger_oring.csv',',',1,0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 134,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(oring(:,2),oring(:,3),'o')\n",
+ "xlabel('Temp (F)')\n",
+ "ylabel('failure (1)/ pass (0)')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 135,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "function J=sse_logistic(a,x,y)\n",
+ " % Create function to calculate cost of logistic function\n",
+ " % t = a0+a1*x\n",
+ " % sigma(t) = 1./(1+e^(-t))\n",
+ " sigma=@(t) 1./(1+exp(-t));\n",
+ " a0=a(1);\n",
+ " a1=a(2);\n",
+ " t=a0+a1*x;\n",
+ " J = 1/length(x)*sum(-y.*log(sigma(t))-(1-y).*log(1-sigma(t)));\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 142,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "J = 0.88822\n",
+ "a =\n",
+ "\n",
+ " 15.03501 -0.23205\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "J=sse_logistic([10,-0.2],oring(:,2),oring(:,3))\n",
+ "a=fminsearch(@(a) sse_logistic(a,oring(:,2),oring(:,3)),[0,-3])\n",
+ "\n",
+ "T=linspace(50,85);\n",
+ "plot(oring(:,2),oring(:,3),'o',T,sigma(a(1)+a(2)*T),T,a(1)+a(2)*T)\n",
+ "axis([50,85,-0.1,1.2])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 139,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "probability of failure when 70 degrees is 23.00% \n",
+ "probability of failure when 60 degrees is 75.25%\n",
+ "probability of failure when 36 degrees is 99.87%\n"
+ ]
+ }
+ ],
+ "source": [
+ "fprintf('probability of failure when 70 degrees is %1.2f%% ',100*sigma(a(1)+a(2)*70))\n",
+ "fprintf('probability of failure when 60 degrees is %1.2f%%',100*sigma(a(1)+a(2)*60))\n",
+ "fprintf('probability of failure when 36 degrees is %1.2f%%',100*sigma(a(1)+a(2)*36))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Interpolation\n",
+ "\n",
+ "Using regression (linear and nonlinear) you are faced with the problem, that you have lots of noisy data and you want to fit a physical model to it. \n",
+ "\n",
+ "You can use interpolation to solve the opposite problem, you have a little data with very little noise.\n",
+ "\n",
+ "## Linear interpolation\n",
+ "\n",
+ "If you are trying to find the value of f(x) for x between $x_{1}$ and $x_{2}$, then you can match the slopes\n",
+ "\n",
+ "$\\frac{f(x)-f(x_{1})}{x-x_{1}}=\\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}$\n",
+ "\n",
+ "or\n",
+ "\n",
+ "$f(x)=f(x_{1})+(x-x_{1})\\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}$\n",
+ "\n",
+ "### Example: Logarithms\n",
+ "\n",
+ "Engineers used to have to use interpolation in logarithm tables for calculations. Find ln(2) from \n",
+ "\n",
+ "a. ln(1) and ln(6)\n",
+ "\n",
+ "b. ln(1) and ln(4)\n",
+ "\n",
+ "c. just calculate it as ln(2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 149,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ln(2)~0.358352\n",
+ "ln(2)~0.462098\n",
+ "ln(2)~0.549306\n",
+ "ln(2)=0.693147\n"
+ ]
+ }
+ ],
+ "source": [
+ "ln2_16=log(1)+(log(6)-log(1))/(6-1)*(2-1);\n",
+ "fprintf('ln(2)~%f\\n',ln2_16)\n",
+ "ln2_14=log(1)+(log(4)-log(1))/(4-1)*(2-1);\n",
+ "ln2_13=log(1)+(log(3)-log(1))/(3-1)*(2-1);\n",
+ "fprintf('ln(2)~%f\\n',ln2_14)\n",
+ "fprintf('ln(2)~%f\\n',ln2_13)\n",
+ "ln2=log(2);\n",
+ "fprintf('ln(2)=%f\\n',ln2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 147,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(1,6);\n",
+ "plot(x,log(x),2,log(2),'*',...\n",
+ "[1,2,6],[log(1),ln2_16,log(6)],'o-',...\n",
+ "[1,2,4],[log(1),ln2_14,log(4)],'s-')\n",
+ "ylabel('ln(x)')\n",
+ "xlabel('x')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Quadratic interpolation (intro curvature)\n",
+ "\n",
+ "Assume function is parabola between 3 points. The function is can be written as:\n",
+ "\n",
+ "$f_{2}(x)=b_{1}+b_{2}(x-x_{1})+b_{3}(x-x_{1})(x-x_{2})$\n",
+ "\n",
+ "When $x=x_{1}$\n",
+ "\n",
+ "$f(x_{1})=b_{1}$\n",
+ "\n",
+ "when $x=x_{2}$\n",
+ "\n",
+ "$b_{2}=\\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}$\n",
+ "\n",
+ "when $x=x_{3}$\n",
+ "\n",
+ "$b_{3}=\\frac{\\frac{f(x_{3})-f(x_{2})}{x_{3}-x_{2}}\n",
+ "-\\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}}{x_{3}-x_{1}}$\n",
+ "\n",
+ "#### Reexamining the ln(2) with ln(1), ln(4), and ln(6):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 154,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Z =\n",
+ "\n",
+ " 1 1 1\n",
+ " 1 4 16\n",
+ " 1 600 360000\n",
+ "\n",
+ "ans = 5.1766e+05\n",
+ "ans =\n",
+ "\n",
+ " -4.6513e-01\n",
+ " 4.6589e-01\n",
+ " -7.5741e-04\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x=[1,4,600]';\n",
+ "Z=[x.^0,x.^1,x.^2]\n",
+ "cond(Z)\n",
+ "Z\\log(x)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 155,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "b1 = 0\n",
+ "b2 = 0.46210\n",
+ "b3 = -0.051873\n"
+ ]
+ }
+ ],
+ "source": [
+ "x1=1;\n",
+ "x2=4;\n",
+ "x3=6;\n",
+ "f1=log(x1);\n",
+ "f2=log(x2);\n",
+ "f3=log(x3);\n",
+ "\n",
+ "b1=f1\n",
+ "b2=(f2-b1)/(x2-x1)\n",
+ "b3=(f3-f2)/(x3-x2)-b2;\n",
+ "b3=b3/(x3-x1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 157,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 0.56584\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(1,6);\n",
+ "f=@(x) b1+b2*(x-x1)+b3*(x-x1).*(x-x2);\n",
+ "plot(x,log(x),2,log(2),'*',...\n",
+ "[1,4,6],[log(1),log(4),log(6)],'ro',...\n",
+ "x,f(x),'r-',2,f(2),'s')\n",
+ "f(2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Newton's Interpolating Polynomials\n",
+ "\n",
+ "For n-data points, we can fit an (n-1)th-polynomial\n",
+ "\n",
+ "$f_{n-1}(x)=b_{1}+b_{2}(x-x_{1})+\\cdots+b_{n}(x-x_{1})(x-x_{2})\\cdots(x-x_{n})$\n",
+ "\n",
+ "where \n",
+ "\n",
+ "$b_{1}=f(x_{1})$\n",
+ "\n",
+ "$b_{2}=\\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}$\n",
+ "\n",
+ "$b_{3}=\\frac{\\frac{f(x_{3})-f(x_{2})}{x_{3}-x_{2}}\n",
+ "-b_{2}}{x_{3}-x_{1}}$\n",
+ "\n",
+ "$\\vdots$\n",
+ "\n",
+ "$b_{n}=f(x_{n},x_{n-1},...,x_{2},x_{1})\n",
+ "=\\frac{f(x_{n},x_{n-1},...x_{2})-f(x_{n-1},x_{n-2},...,x_{1})}{x_{n}-x_{1}}$\n",
+ "\n",
+ "**e.g. for 4 data points:**\n",
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 160,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "b =\n",
+ "\n",
+ " 0.00000 0.54931 -0.08721 0.01178\n",
+ " 1.09861 0.28768 -0.02832 0.00000\n",
+ " 1.38629 0.20273 0.00000 0.00000\n",
+ " 1.79176 0.00000 0.00000 0.00000\n",
+ "\n",
+ "ans = 0.66007\n",
+ "ans =\n",
+ "\n",
+ " 0.00000\n",
+ " 1.09861\n",
+ " 1.38629\n",
+ " 1.79176\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "Newtint([1,3,4,6],log([1,3,4,6]),2)\n",
+ "log([1,3,4,6]')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ln(2)=0.693147\n",
+ "ln(2)~0.693147\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "\n",
+ "\n",
+ "x=[0.2,1,2,3,4]; % define independent var's\n",
+ "y=log(x); % define dependent var's\n",
+ "xx=linspace(min(x),max(x)+1);\n",
+ "yy=zeros(size(xx));\n",
+ "for i=1:length(xx)\n",
+ " yy(i)=Newtint(x,y,xx(i));\n",
+ "end\n",
+ "plot(xx,log(xx),2,log(2),'*',...\n",
+ "x,y,'ro',...\n",
+ "xx,yy,'r-')\n",
+ "\n",
+ "fprintf('ln(2)=%f',log(2))\n",
+ "fprintf('ln(2)~%f',Newtint(x,y,2))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_18/newton_interpolation.png b/lecture_18/newton_interpolation.png
new file mode 100644
index 0000000..5990cb5
Binary files /dev/null and b/lecture_18/newton_interpolation.png differ
diff --git a/lecture_18/octave-workspace b/lecture_18/octave-workspace
new file mode 100644
index 0000000..d40f9ac
Binary files /dev/null and b/lecture_18/octave-workspace differ
diff --git a/lecture_19/.ipynb_checkpoints/lecture 19-checkpoint.ipynb b/lecture_19/.ipynb_checkpoints/lecture 19-checkpoint.ipynb
new file mode 100644
index 0000000..2911efe
--- /dev/null
+++ b/lecture_19/.ipynb_checkpoints/lecture 19-checkpoint.ipynb
@@ -0,0 +1,497 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Splines (Brief introduction before next section)\n",
+ "\n",
+ "Following interpolation discussion, instead of estimating 9 data points with an eighth-order polynomial, it makes more sense to fit sections of the curve to lower-order polynomials:\n",
+ "\n",
+ "0. zeroth-order (nearest neighbor)\n",
+ "\n",
+ "1. first-order (linear interpolation)\n",
+ "\n",
+ "2. third-order (cubic interpolation)\n",
+ "\n",
+ "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n",
+ "\n",
+ "`interp1`\n",
+ "\n",
+ "`interp2`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n",
+ "\n",
+ " -- Function File: YI = interp1 (X, Y, XI)\n",
+ " -- Function File: YI = interp1 (Y, XI)\n",
+ " -- Function File: YI = interp1 (..., METHOD)\n",
+ " -- Function File: YI = interp1 (..., EXTRAP)\n",
+ " -- Function File: YI = interp1 (..., \"left\")\n",
+ " -- Function File: YI = interp1 (..., \"right\")\n",
+ " -- Function File: PP = interp1 (..., \"pp\")\n",
+ "\n",
+ " One-dimensional interpolation.\n",
+ "\n",
+ " Interpolate input data to determine the value of YI at the points\n",
+ " XI. If not specified, X is taken to be the indices of Y ('1:length\n",
+ " (Y)'). If Y is a matrix or an N-dimensional array, the\n",
+ " interpolation is performed on each column of Y.\n",
+ "\n",
+ " The interpolation METHOD is one of:\n",
+ "\n",
+ " \"nearest\"\n",
+ " Return the nearest neighbor.\n",
+ "\n",
+ " \"previous\"\n",
+ " Return the previous neighbor.\n",
+ "\n",
+ " \"next\"\n",
+ " Return the next neighbor.\n",
+ "\n",
+ " \"linear\" (default)\n",
+ " Linear interpolation from nearest neighbors.\n",
+ "\n",
+ " \"pchip\"\n",
+ " Piecewise cubic Hermite interpolating\n",
+ " polynomial--shape-preserving interpolation with smooth first\n",
+ " derivative.\n",
+ "\n",
+ " \"cubic\"\n",
+ " Cubic interpolation (same as \"pchip\").\n",
+ "\n",
+ " \"spline\"\n",
+ " Cubic spline interpolation--smooth first and second\n",
+ " derivatives throughout the curve.\n",
+ "\n",
+ " Adding '*' to the start of any method above forces 'interp1' to\n",
+ " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n",
+ " referenced. This is usually faster, and is never slower. The\n",
+ " default method is \"linear\".\n",
+ "\n",
+ " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n",
+ " the endpoints using the current METHOD. If EXTRAP is a number,\n",
+ " then replace values beyond the endpoints with that number. When\n",
+ " unspecified, EXTRAP defaults to 'NA'.\n",
+ "\n",
+ " If the string argument \"pp\" is specified, then XI should not be\n",
+ " supplied and 'interp1' returns a piecewise polynomial object. This\n",
+ " object can later be used with 'ppval' to evaluate the\n",
+ " interpolation. There is an equivalence, such that 'ppval (interp1\n",
+ " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n",
+ "\n",
+ " Duplicate points in X specify a discontinuous interpolant. There\n",
+ " may be at most 2 consecutive points with the same value. If X is\n",
+ " increasing, the default discontinuous interpolant is\n",
+ " right-continuous. If X is decreasing, the default discontinuous\n",
+ " interpolant is left-continuous. The continuity condition of the\n",
+ " interpolant may be specified by using the options \"left\" or \"right\"\n",
+ " to select a left-continuous or right-continuous interpolant,\n",
+ " respectively. Discontinuous interpolation is only allowed for\n",
+ " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n",
+ " must be unique.\n",
+ "\n",
+ " An example of the use of 'interp1' is\n",
+ "\n",
+ " xf = [0:0.05:10];\n",
+ " yf = sin (2*pi*xf/5);\n",
+ " xp = [0:10];\n",
+ " yp = sin (2*pi*xp/5);\n",
+ " lin = interp1 (xp, yp, xf);\n",
+ " near = interp1 (xp, yp, xf, \"nearest\");\n",
+ " pch = interp1 (xp, yp, xf, \"pchip\");\n",
+ " spl = interp1 (xp, yp, xf, \"spline\");\n",
+ " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n",
+ " xp,yp,\"r*\");\n",
+ " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n",
+ "\n",
+ " See also: pchip, spline, interpft, interp2, interp3, interpn.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help interp1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'interp2' is a function from the file /usr/share/octave/4.0.0/m/general/interp2.m\n",
+ "\n",
+ " -- Function File: ZI = interp2 (X, Y, Z, XI, YI)\n",
+ " -- Function File: ZI = interp2 (Z, XI, YI)\n",
+ " -- Function File: ZI = interp2 (Z, N)\n",
+ " -- Function File: ZI = interp2 (Z)\n",
+ " -- Function File: ZI = interp2 (..., METHOD)\n",
+ " -- Function File: ZI = interp2 (..., METHOD, EXTRAP)\n",
+ "\n",
+ " Two-dimensional interpolation.\n",
+ "\n",
+ " Interpolate reference data X, Y, Z to determine ZI at the\n",
+ " coordinates XI, YI. The reference data X, Y can be matrices, as\n",
+ " returned by 'meshgrid', in which case the sizes of X, Y, and Z must\n",
+ " be equal. If X, Y are vectors describing a grid then 'length (X)\n",
+ " == columns (Z)' and 'length (Y) == rows (Z)'. In either case the\n",
+ " input data must be strictly monotonic.\n",
+ "\n",
+ " If called without X, Y, and just a single reference data matrix Z,\n",
+ " the 2-D region 'X = 1:columns (Z), Y = 1:rows (Z)' is assumed.\n",
+ " This saves memory if the grid is regular and the distance between\n",
+ " points is not important.\n",
+ "\n",
+ " If called with a single reference data matrix Z and a refinement\n",
+ " value N, then perform interpolation over a grid where each original\n",
+ " interval has been recursively subdivided N times. This results in\n",
+ " '2^N-1' additional points for every interval in the original grid.\n",
+ " If N is omitted a value of 1 is used. As an example, the interval\n",
+ " [0,1] with 'N==2' results in a refined interval with points at [0,\n",
+ " 1/4, 1/2, 3/4, 1].\n",
+ "\n",
+ " The interpolation METHOD is one of:\n",
+ "\n",
+ " \"nearest\"\n",
+ " Return the nearest neighbor.\n",
+ "\n",
+ " \"linear\" (default)\n",
+ " Linear interpolation from nearest neighbors.\n",
+ "\n",
+ " \"pchip\"\n",
+ " Piecewise cubic Hermite interpolating\n",
+ " polynomial--shape-preserving interpolation with smooth first\n",
+ " derivative.\n",
+ "\n",
+ " \"cubic\"\n",
+ " Cubic interpolation (same as \"pchip\").\n",
+ "\n",
+ " \"spline\"\n",
+ " Cubic spline interpolation--smooth first and second\n",
+ " derivatives throughout the curve.\n",
+ "\n",
+ " EXTRAP is a scalar number. It replaces values beyond the endpoints\n",
+ " with EXTRAP. Note that if EXTRAPVAL is used, METHOD must be\n",
+ " specified as well. If EXTRAP is omitted and the METHOD is\n",
+ " \"spline\", then the extrapolated values of the \"spline\" are used.\n",
+ " Otherwise the default EXTRAP value for any other METHOD is \"NA\".\n",
+ "\n",
+ " See also: interp1, interp3, interpn, meshgrid.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help interp2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(-pi,pi,9);\n",
+ "xi=linspace(-pi,pi,100);\n",
+ "y=sin(x);\n",
+ "yi_lin=interp1(x,y,xi,'linear');\n",
+ "yi_spline=interp1(x,y,xi,'spline'); \n",
+ "yi_cubic=interp1(x,y,xi,'cubic');\n",
+ "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n",
+ "axis([-pi,pi,-1.5,1.5])\n",
+ "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example: Accelerate then hold velocity\n",
+ "\n",
+ "Here the time is given as vector t in seconds and the "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "t=[0 2 40 56 68 80 84 96 104 110]';\n",
+ "v=[0 20 20 38 80 80 100 100 125 125]';\n"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_19/Newtint.m b/lecture_19/Newtint.m
new file mode 100644
index 0000000..e4c6c83
--- /dev/null
+++ b/lecture_19/Newtint.m
@@ -0,0 +1,34 @@
+function yint = Newtint(x,y,xx)
+% Newtint: Newton interpolating polynomial
+% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton
+% interpolating polynomial based on n data points (x, y)
+% to determine a value of the dependent variable (yint)
+% at a given value of the independent variable, xx.
+% input:
+% x = independent variable
+% y = dependent variable
+% xx = value of independent variable at which
+% interpolation is calculated
+% output:
+% yint = interpolated value of dependent variable
+
+% compute the finite divided differences in the form of a
+% difference table
+n = length(x);
+if length(y)~=n, error('x and y must be same length'); end
+b = zeros(n,n);
+% assign dependent variables to the first column of b.
+b(:,1) = y(:); % the (:) ensures that y is a column vector.
+for j = 2:n
+ for i = 1:n-j+1
+ b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
+ end
+end
+%b
+% use the finite divided differences to interpolate
+xt = 1;
+yint = b(1,1);
+for j = 1:n-1
+ xt = xt*(xx-x(j));
+ yint = yint+b(1,j+1)*xt;
+end
diff --git a/lecture_19/lecture 19.ipynb b/lecture_19/lecture 19.ipynb
new file mode 100644
index 0000000..fde5c6e
--- /dev/null
+++ b/lecture_19/lecture 19.ipynb
@@ -0,0 +1,1842 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 69,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 70,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Questions from last class\n",
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 74,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " 1.04210\n",
+ " 8.19609\n",
+ " 0.50283\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% for linear regression\n",
+ "x=linspace(0,20)';\n",
+ "y=x.^2 +10*exp(-2*x)+0.5*sinh(x/2)+rand(size(x))*200-100;\n",
+ "Z=[x.^2 exp(-2*x) sinh(x/2)];\n",
+ "a=Z\\y\n",
+ "\n",
+ "plot(x,y,'o',x,a(1)*x.^2+a(2)*exp(-2*x)+a(3)*sinh(x/2))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "# Answer: It depends\n",
+ "\n",
+ "#### Other:\n",
+ "\n",
+ "Twice the amount of points needed\n",
+ "\n",
+ "depends on what order polynomial it is and how far the data needs to be extrapolated \n",
+ "\n",
+ "As man you as possible \n",
+ "\n",
+ "Never extrapolate unless linear interpolation.\n",
+ "\n",
+ "You shouldn't. 2 if the linear is a good fit for the region, and you absolutely have to.\n",
+ "\n",
+ "Wait can you do that?\n",
+ "\n",
+ "Don't use extrapolation\n",
+ "\n",
+ "do not extrapolate\n",
+ "\n",
+ "As many data points as you have\n",
+ "\n",
+ "the more the better so that the best polynomial can be made through the data points\n",
+ "\n",
+ "Twice the amount of points needed"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Questions from you\n",
+ "\n",
+ "- when will the project assignment be finalized? Also do you pronounce it \"jiff\" or \"gif\"?\n",
+ "\n",
+ "- If blue is red and red is blue, then what is purple? \n",
+ "\n",
+ "- How do we open the .ipynb lecture files? Or will the lectures continue to be also saved in pdf (last few have not).\n",
+ "\n",
+ "- When will we be put on teams for the final project?\n",
+ "\n",
+ "- What is the grading rubric for the project?\n",
+ "\n",
+ "- How to sync repository with files from laptop like hw without using Github desktop \n",
+ "\n",
+ "- Are there any upcoming deadlines for the project?\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Splines (Brief introduction before integrals)\n",
+ "\n",
+ "Following interpolation discussion, instead of estimating 9 data points with an eighth-order polynomial, it makes more sense to fit sections of the curve to lower-order polynomials:\n",
+ "\n",
+ "0. zeroth-order (nearest neighbor)\n",
+ "\n",
+ "1. first-order (linear interpolation)\n",
+ "\n",
+ "2. third-order (cubic interpolation)\n",
+ "\n",
+ "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n",
+ "\n",
+ "`interp1`\n",
+ "\n",
+ "`interp2`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n",
+ "\n",
+ " -- Function File: YI = interp1 (X, Y, XI)\n",
+ " -- Function File: YI = interp1 (Y, XI)\n",
+ " -- Function File: YI = interp1 (..., METHOD)\n",
+ " -- Function File: YI = interp1 (..., EXTRAP)\n",
+ " -- Function File: YI = interp1 (..., \"left\")\n",
+ " -- Function File: YI = interp1 (..., \"right\")\n",
+ " -- Function File: PP = interp1 (..., \"pp\")\n",
+ "\n",
+ " One-dimensional interpolation.\n",
+ "\n",
+ " Interpolate input data to determine the value of YI at the points\n",
+ " XI. If not specified, X is taken to be the indices of Y ('1:length\n",
+ " (Y)'). If Y is a matrix or an N-dimensional array, the\n",
+ " interpolation is performed on each column of Y.\n",
+ "\n",
+ " The interpolation METHOD is one of:\n",
+ "\n",
+ " \"nearest\"\n",
+ " Return the nearest neighbor.\n",
+ "\n",
+ " \"previous\"\n",
+ " Return the previous neighbor.\n",
+ "\n",
+ " \"next\"\n",
+ " Return the next neighbor.\n",
+ "\n",
+ " \"linear\" (default)\n",
+ " Linear interpolation from nearest neighbors.\n",
+ "\n",
+ " \"pchip\"\n",
+ " Piecewise cubic Hermite interpolating\n",
+ " polynomial--shape-preserving interpolation with smooth first\n",
+ " derivative.\n",
+ "\n",
+ " \"cubic\"\n",
+ " Cubic interpolation (same as \"pchip\").\n",
+ "\n",
+ " \"spline\"\n",
+ " Cubic spline interpolation--smooth first and second\n",
+ " derivatives throughout the curve.\n",
+ "\n",
+ " Adding '*' to the start of any method above forces 'interp1' to\n",
+ " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n",
+ " referenced. This is usually faster, and is never slower. The\n",
+ " default method is \"linear\".\n",
+ "\n",
+ " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n",
+ " the endpoints using the current METHOD. If EXTRAP is a number,\n",
+ " then replace values beyond the endpoints with that number. When\n",
+ " unspecified, EXTRAP defaults to 'NA'.\n",
+ "\n",
+ " If the string argument \"pp\" is specified, then XI should not be\n",
+ " supplied and 'interp1' returns a piecewise polynomial object. This\n",
+ " object can later be used with 'ppval' to evaluate the\n",
+ " interpolation. There is an equivalence, such that 'ppval (interp1\n",
+ " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n",
+ "\n",
+ " Duplicate points in X specify a discontinuous interpolant. There\n",
+ " may be at most 2 consecutive points with the same value. If X is\n",
+ " increasing, the default discontinuous interpolant is\n",
+ " right-continuous. If X is decreasing, the default discontinuous\n",
+ " interpolant is left-continuous. The continuity condition of the\n",
+ " interpolant may be specified by using the options \"left\" or \"right\"\n",
+ " to select a left-continuous or right-continuous interpolant,\n",
+ " respectively. Discontinuous interpolation is only allowed for\n",
+ " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n",
+ " must be unique.\n",
+ "\n",
+ " An example of the use of 'interp1' is\n",
+ "\n",
+ " xf = [0:0.05:10];\n",
+ " yf = sin (2*pi*xf/5);\n",
+ " xp = [0:10];\n",
+ " yp = sin (2*pi*xp/5);\n",
+ " lin = interp1 (xp, yp, xf);\n",
+ " near = interp1 (xp, yp, xf, \"nearest\");\n",
+ " pch = interp1 (xp, yp, xf, \"pchip\");\n",
+ " spl = interp1 (xp, yp, xf, \"spline\");\n",
+ " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n",
+ " xp,yp,\"r*\");\n",
+ " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n",
+ "\n",
+ " See also: pchip, spline, interpft, interp2, interp3, interpn.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help interp1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'interp2' is a function from the file /usr/share/octave/4.0.0/m/general/interp2.m\n",
+ "\n",
+ " -- Function File: ZI = interp2 (X, Y, Z, XI, YI)\n",
+ " -- Function File: ZI = interp2 (Z, XI, YI)\n",
+ " -- Function File: ZI = interp2 (Z, N)\n",
+ " -- Function File: ZI = interp2 (Z)\n",
+ " -- Function File: ZI = interp2 (..., METHOD)\n",
+ " -- Function File: ZI = interp2 (..., METHOD, EXTRAP)\n",
+ "\n",
+ " Two-dimensional interpolation.\n",
+ "\n",
+ " Interpolate reference data X, Y, Z to determine ZI at the\n",
+ " coordinates XI, YI. The reference data X, Y can be matrices, as\n",
+ " returned by 'meshgrid', in which case the sizes of X, Y, and Z must\n",
+ " be equal. If X, Y are vectors describing a grid then 'length (X)\n",
+ " == columns (Z)' and 'length (Y) == rows (Z)'. In either case the\n",
+ " input data must be strictly monotonic.\n",
+ "\n",
+ " If called without X, Y, and just a single reference data matrix Z,\n",
+ " the 2-D region 'X = 1:columns (Z), Y = 1:rows (Z)' is assumed.\n",
+ " This saves memory if the grid is regular and the distance between\n",
+ " points is not important.\n",
+ "\n",
+ " If called with a single reference data matrix Z and a refinement\n",
+ " value N, then perform interpolation over a grid where each original\n",
+ " interval has been recursively subdivided N times. This results in\n",
+ " '2^N-1' additional points for every interval in the original grid.\n",
+ " If N is omitted a value of 1 is used. As an example, the interval\n",
+ " [0,1] with 'N==2' results in a refined interval with points at [0,\n",
+ " 1/4, 1/2, 3/4, 1].\n",
+ "\n",
+ " The interpolation METHOD is one of:\n",
+ "\n",
+ " \"nearest\"\n",
+ " Return the nearest neighbor.\n",
+ "\n",
+ " \"linear\" (default)\n",
+ " Linear interpolation from nearest neighbors.\n",
+ "\n",
+ " \"pchip\"\n",
+ " Piecewise cubic Hermite interpolating\n",
+ " polynomial--shape-preserving interpolation with smooth first\n",
+ " derivative.\n",
+ "\n",
+ " \"cubic\"\n",
+ " Cubic interpolation (same as \"pchip\").\n",
+ "\n",
+ " \"spline\"\n",
+ " Cubic spline interpolation--smooth first and second\n",
+ " derivatives throughout the curve.\n",
+ "\n",
+ " EXTRAP is a scalar number. It replaces values beyond the endpoints\n",
+ " with EXTRAP. Note that if EXTRAPVAL is used, METHOD must be\n",
+ " specified as well. If EXTRAP is omitted and the METHOD is\n",
+ " \"spline\", then the extrapolated values of the \"spline\" are used.\n",
+ " Otherwise the default EXTRAP value for any other METHOD is \"NA\".\n",
+ "\n",
+ " See also: interp1, interp3, interpn, meshgrid.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help interp2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 75,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(-pi,pi,9);\n",
+ "xi=linspace(-pi,pi,100);\n",
+ "y=sin(x);\n",
+ "yi_lin=interp1(x,y,xi,'linear');\n",
+ "yi_spline=interp1(x,y,xi,'spline'); \n",
+ "yi_cubic=interp1(x,y,xi,'cubic');\n",
+ "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n",
+ "axis([-pi,pi,-1.5,1.5])\n",
+ "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example: Accelerate then hold velocity\n",
+ "\n",
+ "Test driving a car, the accelerator is pressed, then released, then pressed again for 20-second intervals, until speed is 120 mph. Here the time is given as vector t in seconds and the velocity is in mph. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 82,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "t=[0 20 40 68 80 84 96 104 110]';\n",
+ "v=[0 20 20 80 80 100 100 125 125]';\n",
+ "tt=linspace(0,110)';\n",
+ "v_lin=interp1(t,v,tt);\n",
+ "v_spl=interp1(t,v,tt,'spline');\n",
+ "v_cub=interp1(t,v,tt,'cubic');\n",
+ "\n",
+ "plot(t,v,'o',56,38,'s',tt,v_lin,tt,v_spl,tt,v_cub)\n",
+ "xlabel('t (s)')\n",
+ "ylabel('v (mph)')\n",
+ "legend('data','removed data point','linear','cubic spline','piecewise cubic','Location','NorthWest')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 83,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "t=[0 20 40 56 68 80 84 96 104 110]';\n",
+ "v=[0 20 20 38 80 80 100 100 125 125]';\n",
+ "tt=linspace(0,110)';\n",
+ "v_lin=interp1(t,v,tt);\n",
+ "v_spl=interp1(t,v,tt,'spline');\n",
+ "v_cub=interp1(t,v,tt,'cubic');\n",
+ "\n",
+ "\n",
+ "plot(tt(2:end),diff(v_lin)./diff(tt),tt(2:end),diff(v_spl)./diff(tt),tt(2:end),diff(v_cub)./diff(tt))\n",
+ "xlabel('t (s)')\n",
+ "ylabel('dv/dt (mph/s)')\n",
+ "legend('linear','cubic spline','piecewise cubic','Location','NorthWest')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Choose spline wisely\n",
+ "\n",
+ "For example of sin(x), not very important\n",
+ "\n",
+ "For stop-and-hold examples, the $C^{2}$-continuity should not be preserved. You don't need smooth curves.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Numerical Integration\n",
+ "\n",
+ "A definite integral is defined by \n",
+ "\n",
+ "$I=\\int_{a}^{b}f(x)dx$\n",
+ "\n",
+ "To determine the mass of an object with varying density, you can perform a summation\n",
+ "\n",
+ "mass=$\\sum_{i=1}^{n}\\rho_{i}\\Delta V_{i}$\n",
+ "\n",
+ "or taking the limit as $\\Delta V \\rightarrow dV=dxdydz$\n",
+ "\n",
+ "mass=$\\int_{0}^{h}\\int_{0}^{w}\\int_{0}^{l}\\rho(x,y,z)dxdydz$\n",
+ "\n",
+ "## Newton-Cotes Formulas\n",
+ "\n",
+ "$I=\\int_{a}^{b}f(x)dx=\\int_{a}^{b}f_{n}(x)dx$\n",
+ "\n",
+ "where $f_{n}$ is an n$^{th}$-order polynomial approximation of f(x)\n",
+ "\n",
+ "## First-Order: Trapezoidal Rule\n",
+ "\n",
+ "$I=\\int_{a}^{b}f(x)dx\\approx \\int_{a}^{b}\\left(f(a)+\\frac{f(b)-f(a)}{b-a}(x-a)\\right)dx$\n",
+ "\n",
+ "$I\\approx(b-a)\\frac{f(a)+f(b)}{2}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "I_trap = 0.78540\n",
+ "I_act = 1.00000\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(0,pi)';\n",
+ "plot(x,sin(x),[0,pi/2],sin([0,pi/2]))\n",
+ "I_trap=mean(sin(([0,pi/2]))*(diff([0,pi/2])))\n",
+ "I_act = -(cos(pi/2)-cos(0))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Improve estimate with more points\n",
+ "\n",
+ "$I=\\int_{a}^{b}f(x)dx=\\int_{a}^{a+\\Delta x}f(x)dx+\\int_{a+\\Delta x}^{a+2\\Delta x}f(x)dx+ \\cdots \\int_{b-\\Delta x}^{b}f(x)dx$\n",
+ "\n",
+ "$I\\approx\\Delta x\\frac{f(a)+f(a+\\Delta x)}{2}+\\Delta x\\frac{f(a+\\Delta x)+f(a+2\\Delta x)}{2}\n",
+ "+\\cdots \\Delta x\\frac{f(b-\\Delta x)+f(b)}{2}$\n",
+ "\n",
+ "$I\\approx \\frac{\\Delta x}{2}\\left(f(a)+2\\sum_{i=1}^{n-1}f(a+i\\Delta x) +f(b)\\right)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 86,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "For 70 steps\n",
+ "trapezoid approximation of integral is 0.99 \n",
+ " actual integral is 1.00\n"
+ ]
+ }
+ ],
+ "source": [
+ "N=70;\n",
+ "I_trap=trap(@(x) sin(x),0,pi/2,N);\n",
+ "fprintf('For %i steps\\ntrapezoid approximation of integral is %1.2f \\n actual integral is %1.2f',N,I_trap,I_act)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 92,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(0,pi);\n",
+ "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'-o')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Increase accuracy = Increase polynomial order\n",
+ "\n",
+ "### Simpson's Rules\n",
+ "\n",
+ "When integrating f(x) and using a second order polynomial, this is known as **Simpson's 1/3 Rule**\n",
+ "\n",
+ "$I=\\frac{h}{3}(f(x_{0})+4f(x_{1})+f(x_{2}))$\n",
+ "\n",
+ "where a=$x_{0}$, b=$x_{2}$, and $x_{1}=\\frac{a+b}{2}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This can be used with n=3 or multiples of 2 intervals\n",
+ "\n",
+ "$I=\\int_{x_{0}}^{x_{2}}f(x)dx+\\int_{x_{2}}^{x_{4}}f(x)dx+\\cdots +\\int_{x_{n-2}}^{x_{n}}f(x)dx$\n",
+ "\n",
+ "$I=(b-a)\\frac{f(x_{0})+4\\sum_{i=1,3,5}^{n-1}f(x_{i})+2\\sum_{i=2,4,6}^{n-2}f(x_{i})+f(x_{n})}{3n}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 93,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1.6235\n",
+ "Is_1_3 = 1.0023\n"
+ ]
+ }
+ ],
+ "source": [
+ "f=@(x) 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;\n",
+ "simpson3(f,0,0.8,4)\n",
+ "Is_1_3=simpson3(@(x) sin(x),0,pi/2,2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## General Newton-Cotes formulae\n",
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_19/newton_cotes.png b/lecture_19/newton_cotes.png
new file mode 100644
index 0000000..b734d11
Binary files /dev/null and b/lecture_19/newton_cotes.png differ
diff --git a/lecture_19/octave-workspace b/lecture_19/octave-workspace
new file mode 100644
index 0000000..4bd89dd
Binary files /dev/null and b/lecture_19/octave-workspace differ
diff --git a/lecture_19/q1.png b/lecture_19/q1.png
new file mode 100644
index 0000000..44ec9da
Binary files /dev/null and b/lecture_19/q1.png differ
diff --git a/lecture_19/q2.png b/lecture_19/q2.png
new file mode 100644
index 0000000..fb27a54
Binary files /dev/null and b/lecture_19/q2.png differ
diff --git a/lecture_19/q3.png b/lecture_19/q3.png
new file mode 100644
index 0000000..62c4402
Binary files /dev/null and b/lecture_19/q3.png differ
diff --git a/lecture_19/q4.png b/lecture_19/q4.png
new file mode 100644
index 0000000..5457483
Binary files /dev/null and b/lecture_19/q4.png differ
diff --git a/lecture_19/simpson3.m b/lecture_19/simpson3.m
new file mode 100644
index 0000000..2ae0c81
--- /dev/null
+++ b/lecture_19/simpson3.m
@@ -0,0 +1,23 @@
+function I = simpson3(func,a,b,n,varargin)
+% simpson3: composite simpson's 1/3 rule
+% I = simpson3(func,a,b,n,pl,p2,...):
+% composite trapezoidal rule
+% input:
+% func = name of function to be integrated
+% a, b = integration limits
+% n = number of segments (default = 100)
+% pl,p2,... = additional parameters used by func
+% output:
+% I = integral estimate
+if nargin<3,error('at least 3 input arguments required'),end
+if ~(b>a),error('upper bound must be greater than lower'),end
+if nargin<4|isempty(n),n=100;end
+x = a; h = (b - a)/n;
+
+xvals=linspace(a,b,n+1);
+fvals=func(xvals,varargin{:});
+s=fvals(1);
+s = s + 4*sum(fvals(2:2:end-1));
+s = s + 2*sum(fvals(3:2:end-2));
+s = s + fvals(end);
+I = (b - a) * s/(3*n);
diff --git a/lecture_19/trap.m b/lecture_19/trap.m
new file mode 100644
index 0000000..85b8685
--- /dev/null
+++ b/lecture_19/trap.m
@@ -0,0 +1,22 @@
+function I = trap(func,a,b,n,varargin)
+% trap: composite trapezoidal rule quadrature
+% I = trap(func,a,b,n,pl,p2,...):
+% composite trapezoidal rule
+% input:
+% func = name of function to be integrated
+% a, b = integration limits
+% n = number of segments (default = 100)
+% pl,p2,... = additional parameters used by func
+% output:
+% I = integral estimate
+if nargin<3,error('at least 3 input arguments required'),end
+if ~(b>a),error('upper bound must be greater than lower'),end
+if nargin<4|isempty(n),n=100;end
+
+x = a; h = (b - a)/n;
+xvals=linspace(a,b,n);
+fvals=func(xvals,varargin{:});
+s=func(a,varargin{:});
+s = s + 2*sum(fvals(2:n-1));
+s = s + func(b,varargin{:});
+I = (b - a) * s/(2*n);