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_18/lecture_18.ipynb b/lecture_18/lecture_18.ipynb index bf7731c..f1fd1b3 100644 --- a/lecture_18/lecture_18.ipynb +++ b/lecture_18/lecture_18.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 146, + "execution_count": 2, "metadata": { "collapsed": true }, @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 107, + "execution_count": 3, "metadata": { "collapsed": true }, @@ -1959,7 +1959,7 @@ }, { "cell_type": "code", - "execution_count": 162, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -1969,7 +1969,7 @@ "output_type": "stream", "text": [ "ln(2)=0.693147\n", - "ln(2)~0.366349\n" + "ln(2)~0.693147\n" ] }, { @@ -2014,113 +2014,117 @@ "\n", "\n", "\t\n", - "\t\t\n", + "\t\t\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t-2\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0\n", + "\t\t\n", + "\t\t-1.5\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t2\n", + "\t\t\n", + "\t\t-1\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t4\n", + "\t\t\n", + "\t\t-0.5\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t6\n", + "\t\t\n", + "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t8\n", + "\t\t\n", + "\t\t0.5\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t10\n", + "\t\t\n", + "\t\t1\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0\n", + "\t\t\n", + "\t\t1.5\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t10\n", + "\t\t\n", + "\t\t2\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t20\n", + "\t\t\n", + "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t30\n", + "\t\t\n", + "\t\t1\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t40\n", + "\t\t\n", + "\t\t2\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t50\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", "\t\n", "\n", "\n", "\t\t\n", - "\t\t60\n", + "\t\t5\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\n", + "\t\n", "\n", "\n", "\tgnuplot_plot_1a\n", "\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -2144,9 +2148,9 @@ "source": [ "\n", "\n", - "x=[0.2,3,10,20,50,60]; % define independent var's\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));\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", diff --git a/lecture_18/octave-workspace b/lecture_18/octave-workspace index c651696..d40f9ac 100644 Binary files a/lecture_18/octave-workspace and b/lecture_18/octave-workspace differ diff --git a/lecture_19/lecture 19.ipynb b/lecture_19/lecture 19.ipynb index e6f9f76..fde5c6e 100644 --- a/lecture_19/lecture 19.ipynb +++ b/lecture_19/lecture 19.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 69, "metadata": { "collapsed": true }, @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 70, "metadata": { "collapsed": true }, @@ -38,6 +38,283 @@ "![q2](q2.png)" ] }, + { + "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": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t \n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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": {}, @@ -51,6 +328,8 @@ "source": [ "![q4](q4.png)\n", "\n", + "# Answer: It depends\n", + "\n", "#### Other:\n", "\n", "Twice the amount of points needed\n", @@ -101,7 +380,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Splines (Brief introduction before next section)\n", + "# 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", @@ -315,7 +594,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 75, "metadata": { "collapsed": false }, @@ -532,12 +811,12 @@ "source": [ "### Example: Accelerate then hold velocity\n", "\n", - "Here the time is given as vector t in seconds and the velocity is in mph. " + "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": 10, + "execution_count": 82, "metadata": { "collapsed": false }, @@ -683,12 +962,12 @@ "\n", "\n", "\n", - "\t\n", + "\t\n", "\tdata\n", "\n", "\n", "\n", - "\t\n", + "\t\n", "\t\tdata\n", "\t\n", "\n", @@ -696,7 +975,6 @@ "\t\n", "\t\n", "\t\n", - "\t\n", "\t\n", "\t\n", "\t\n", @@ -706,34 +984,45 @@ "\t\n", "\n", "\t\n", - "\tlinear\n", + "\tremoved data point\n", "\n", - "\t\n", + "\t\n", + "\t\tremoved data point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", "\t\tlinear\n", "\t\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", - "\tcubic spline\n", + "\tcubic spline\n", "\n", - "\t\n", + "\t\n", "\t\tcubic spline\n", "\t\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", - "\tpiecewise cubic\n", + "\tpiecewise cubic\n", "\n", - "\t\n", + "\t\n", "\t\tpiecewise cubic\n", "\t\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", - "\n", + "\n", "\n", "\n", "\n", @@ -753,22 +1042,22 @@ } ], "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", + "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',tt,v_lin,tt,v_spl,tt,v_cub)\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','linear','cubic spline','piecewise cubic','Location','NorthWest')" + "legend('data','removed data point','linear','cubic spline','piecewise cubic','Location','NorthWest')" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 83, "metadata": { "collapsed": false }, @@ -1206,7 +1495,7 @@ }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 86, "metadata": { "collapsed": false }, @@ -1215,21 +1504,21 @@ "name": "stdout", "output_type": "stream", "text": [ - "For 5 steps\n", - "trapezoid approximation of integral is 0.79 \n", + "For 70 steps\n", + "trapezoid approximation of integral is 0.99 \n", " actual integral is 1.00\n" ] } ], "source": [ - "N=5;\n", + "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": 43, + "execution_count": 92, "metadata": { "collapsed": false }, @@ -1365,11 +1654,76 @@ "\t\n", "\tgnuplot_plot_2a\n", "\n", - "\t \n", + "\t\t \n", "\t\n", - "\t\n", - "\t\n", - "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", "\t\n", "\n", "\t\n", @@ -1393,8 +1747,8 @@ } ], "source": [ - "\n", - "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'o')" + "x=linspace(0,pi);\n", + "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'-o')" ] }, { @@ -1425,7 +1779,7 @@ }, { "cell_type": "code", - "execution_count": 68, + "execution_count": 93, "metadata": { "collapsed": false }, 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