Update

HW6/.~lock.Primary_Energy_monthly.csv#

@@ -0,0 +1 @@
+,ryan,fermi,31.03.2017 16:47,file:///home/ryan/.config/libreoffice/4;

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. 
+
+
+

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

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
+

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];
+

lecture_16/lecture_16.ipynb

@@ -20644,7 +20644,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 123,
+   "execution_count": 200,
    "metadata": {
     "collapsed": false
    },
@@ -20793,24 +20793,38 @@
        "</g>\n",
        "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
        "</g>\n",
-       "\t<g id=\"gnuplot_plot_1a\"><title>gnuplot_plot_1a</title>\n",
-       "<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
+       "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
+       "</g>\n",
+       "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
+       "\t<path d=\"M80.2,121.7 L80.2,25.7 L246.0,25.7 L246.0,121.7 L80.2,121.7 Z  \" stroke=\"black\"/></g>\n",
+       "\t<g id=\"gnuplot_plot_1a\"><title>data</title>\n",
+       "<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
        "</g>\n",
        "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
-       "\t<g onmousemove=\"gnuplot_svg.showHypertext(evt,'')\" onmouseout=\"gnuplot_svg.hideHypertext()\"><title> </title>\n",
-       "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(71.9,289.8) scale(12.00)\" xlink:href=\"#gpPt5\"/></g>\n",
+       "\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(246.0,55.7)\">\n",
+       "\t\t<text><tspan font-family=\"{}\">data</tspan></text>\n",
+       "\t</g>\n",
+       "</g>\n",
+       "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
+       "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(71.9,289.8) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(138.1,283.6) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(204.2,240.7) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(270.4,217.2) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(336.5,208.9) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(402.7,124.6) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(468.8,178.5) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(535.0,92.8) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
+       "\t<use color=\"rgb(  0,   0, 255)\" transform=\"translate(118.3,49.7) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
        "</g>\n",
        "\t</g>\n",
-       "\t<g id=\"gnuplot_plot_2a\"><title>gnuplot_plot_2a</title>\n",
+       "\t<g id=\"gnuplot_plot_2a\"><title>best-fit</title>\n",
+       "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
+       "\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(246.0,103.7)\">\n",
+       "\t\t<text><tspan font-family=\"{}\">best-fit</tspan></text>\n",
+       "\t</g>\n",
+       "</g>\n",
        "<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
-       "\t<path d=\"M71.9,298.7 L138.1,271.8 L204.2,244.9 L270.4,218.0 L336.5,191.0 L402.7,164.1 L468.8,137.2 L535.0,110.3    \" stroke=\"rgb(  0, 128,   0)\"/></g>\n",
+       "\t<path d=\"M91.4,97.7 L145.2,97.7 M71.9,298.7 L138.1,271.8 L204.2,244.9 L270.4,218.0 L336.5,191.0 L402.7,164.1   L468.8,137.2 L535.0,110.3  \" stroke=\"rgb(  0, 128,   0)\"/></g>\n",
        "\t</g>\n",
        "<g color=\"white\" fill=\"none\" stroke=\"rgb(  0, 128,   0)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
        "</g>\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)')"
    ]