From 3b1fa91d933d842683684373f06113416242d810 Mon Sep 17 00:00:00 2001 From: ryan Date: Tue, 15 Oct 2019 17:51:13 +0000 Subject: [PATCH] added challenger case study to 15 --- 15_challenger-logistic/.Newtint.m.swp | Bin 0 -> 12288 bytes ...tic-regression_challenger-checkpoint.ipynb | 6 + .../15_logistic-regression_challenger.ipynb | 2197 +++++++++++++++++ 15_challenger-logistic/Newtint.m | 34 + 15_challenger-logistic/challenger_oring.csv | 24 + .../newton_interpolation.png | Bin 0 -> 16742 bytes 15_challenger-logistic/octave-workspace | Bin 0 -> 153 bytes 7 files changed, 2261 insertions(+) create mode 100644 15_challenger-logistic/.Newtint.m.swp create mode 100644 15_challenger-logistic/.ipynb_checkpoints/15_logistic-regression_challenger-checkpoint.ipynb create mode 100644 15_challenger-logistic/15_logistic-regression_challenger.ipynb create mode 100644 15_challenger-logistic/Newtint.m create mode 100644 15_challenger-logistic/challenger_oring.csv create mode 100644 15_challenger-logistic/newton_interpolation.png create mode 100644 15_challenger-logistic/octave-workspace diff --git a/15_challenger-logistic/.Newtint.m.swp b/15_challenger-logistic/.Newtint.m.swp new file mode 100644 index 0000000000000000000000000000000000000000..769fe2bbba029fa7db9ea0d02c35e1f0f2bec202 GIT binary patch literal 12288 zcmeI2KW`K{7{(nC6@LyL(p3~Mpokf=o5>^~U?W5UNKkM>Lb;^?pGZe;g)EuG_pRm2q{lcFh`-+`h50 zeBnY;s?e<8U0zBrCb?~@TCESXOGYYpO{Lb;OUcjbq0`1Kmj7Ta4q*aJV1_{B`rS*% z#L9BI6nf8`K1C;fT$>5tS0=y&m;e)C0!)AjFaajO1pe;?+~%Bkf+7#Jg`R5Lxo_I` zD>qDl2`~XBzyz286JP>NfC(@GCcp%kz<-c{EQI*=un>Ph;s5{JfB%2{L5SDTBj_r0 z2zq}^h*q7tPr(;IOn?b60Vco%m;e)C0!)Aj{Cfh*C{qZcF>ce+l}WRZMQLYaA&t*@ zw~8Sj2E02LdEfQXMjWq-zO85jy=l*MQ_<#zv=oDr_;gjwQ&TH)Lq&aUv{O{*EnTPr zZ6Ei9mC032&JtEu)utUur^M9E^phxybY~;RIuYI$V5*zO-AN~2;w19uo%G=)*dhwX zbZZPH>=<3{8IeUj_`~2h>cy)xKe2{aQl`eWVt5%!M`O}8Ny_c0DGhDm%CW5$QH!i= zJupE1B1FxBV`6|7>(n{S4@Xg zRoN;!>q(jdOP391$Q0GGRBhy0SPeW3rfPFNW}2U&NG|w5S8mJA2B7z0uFW*URwcow z6iD4U0W-T%pvF<7Ld@Q>kz*T&<$9BBOh7X`sjUedb2-YJu(8D?$a@sAGp7c}Z($O4 z;Nei`!>J{N?~R#Gg`w|FfI)3`Twyjz1HGlNliQ^PIN!_ckoRGuDV)hd;Y^m=DCF-=HKwN84%9&*w|9b5oNN*!3PwaAn7EA5L#2@J2Vcd?N6w^){gFR$NRWP wYcHq)$ literal 0 HcmV?d00001 diff --git a/15_challenger-logistic/.ipynb_checkpoints/15_logistic-regression_challenger-checkpoint.ipynb b/15_challenger-logistic/.ipynb_checkpoints/15_logistic-regression_challenger-checkpoint.ipynb new file mode 100644 index 0000000..2fd6442 --- /dev/null +++ b/15_challenger-logistic/.ipynb_checkpoints/15_logistic-regression_challenger-checkpoint.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/15_challenger-logistic/15_logistic-regression_challenger.ipynb b/15_challenger-logistic/15_logistic-regression_challenger.ipynb new file mode 100644 index 0000000..f1fd1b3 --- /dev/null +++ b/15_challenger-logistic/15_logistic-regression_challenger.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": [ + "\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-50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\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", + "\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": [ + "[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": [ + "\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-50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t150\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t200\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t250\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t300\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\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", + "\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": [ + "[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": [ + "\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-15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tresiduals of function\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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": [ + "\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\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", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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": [ + "\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t55\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t65\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t70\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t75\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t85\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tfailure (1)/ pass (0)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tTemp (F)\n", + "\t\n", + "\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", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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": [ + "\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t50\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t55\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t65\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t70\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t75\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t85\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", + "\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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": [ + "\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\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\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tln(x)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tx\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t\t \n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t\t \n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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": [ + "\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\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\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_5a\n", + "\n", + "\t \n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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", + "![Newton Interpolation Iterations](newton_interpolation.png)" + ] + }, + { + "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": [ + "\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-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\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\t5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\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", + "\n", + "\t\n", + "\tgnuplot_plot_4a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "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/15_challenger-logistic/Newtint.m b/15_challenger-logistic/Newtint.m new file mode 100644 index 0000000..e4c6c83 --- /dev/null +++ b/15_challenger-logistic/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/15_challenger-logistic/challenger_oring.csv b/15_challenger-logistic/challenger_oring.csv new file mode 100644 index 0000000..11d647e --- /dev/null +++ b/15_challenger-logistic/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/15_challenger-logistic/newton_interpolation.png b/15_challenger-logistic/newton_interpolation.png new file mode 100644 index 0000000000000000000000000000000000000000..5990cb5fcb1a38d385462ce02d50e804d9f0b2fe GIT binary patch literal 16742 zcmcJ%1yq#nyElq~sFc!;DAL`H(%m5)Lw6&Mf^;cTQqs~O9RngTbjQ%$Dc$hh<9q%) z&i?jUXMgJ)*AfAlXP&t3>;Bb6sFH#t`ctB(NJvQN(o$k7NJ#hHz~>c@AAye#_mts? zf1E|6)gFUC-j7W}z}JK>;_qBk?af@=jhsx8%#nPQa`)>14#lrcpcQpofVC@CFU;|euyIbe2F3gMX_#5Y_t~4qqtom7 zRv2l-o0%qYp`+*oh3Tj%J&*DE&~tP1Vdf!|Lrb_W)BgJMqG;Xv%)rfz^W53tDLbww zR6qKe$YY8?i-WIydiN*-!52nkUnZU)K702>dI#|hilRVDH^dJ;(-6Vp|NAraMv!i4 z^D0j#S<}NJ#E*zMJ5{LZrwy6*F88>IVLV)wcS`=nsiX=_w%;1~G70s1`!jX(3Quao z3yGiN==OHq|8$=_x|*|}mY&ge)IGDRBZJZnZ{s)=l1;VO+ONfMD{edm|B-;bP{mil=6x1{@wr0OaLq=A=6i4Z-mx!N`&N87Gw@BzSf={9ax&JZ zM>-Ti*0lL(ctl~;^Llo_+PF}I-TlyF^c2GzKknUfN{+Wrai=h~Hl+@AjCXfgjy>u(Ro#Rmd4YktE{@cTDsSpyD^ zg_ou|^z_S6d7ZZzZ!EJi#R_v?Fyr$_K9Zz1XEHNacQK&BN#XVQ=xnYqlCR?~!0#qf z3XgTBh)tcY`I6@L@%Kl8pX>w60Y1tVd&K0 z`Bg$BH^7ZUdacbUiJ3bH>r2H*Z5TgS=p5)-ooKnKJlpzc-l3o#6_>Cn#qS&Vu~FCa z_@rg7W$4d)Q*JF{Zuj9RG8k+euhdSkK%iNyaJBP;37^Z}lv;vFGlcH_tK%r-CBuCUx6m4`ySeG>uP!SEP0U6@rlvVTe)i`J zq%b)shCR7((`2mDYGa^a!~28J$SbQR4+iVJ6w^Gt&{G@E9vV*P+P~l4e7(UG>Hdx5 zNP^6J%LC^ra&RoRhvkqTTC?HqsMVU#iyXksO$A&{+BJGQIZ0#7ZgP9bJ>i0}81Z2R#ENmsig9$F#Fv=UFw?cP5F89QdF9gbVv`U&G zOnFh{t0&atC<{ADs*QEtPu#4~&lvTCB(E@RZZ2hSCvQ>P#l;1h|J>_O5okG{-dp;R zmnK?VzMhI7AB$4~SBf6~f|W8=x%k0!AhW#2d?3BK;;LpBQIqqMWZ5r-qqUsZ9nXHZ zlT6^dpD=&5IFON*ad2_We-0TVyi#E39I$nDv`}Vr+uIeWb6$kl%wID$+f7M1IzG4e zt3U=@TI!RVUB)<(48s-9k&v14Nb9wk?g;+%`wlNCFrfO1#N3$Q7vUXavjL!SZFF}l zO}-l~#yi;<$|Eo6sh^8<&C`sM9ojx=Jla)s*`5vVOZU0|D=D^v*Uds$vGJ0s=`Aj{ zrJdJ?5x72usiD?D3{uMw_8~UfIGK3pEt9g}T}!lqkAm)sLVub9 z)ZuRqSKzx}&u1Z#C`d1zkvH3?p;BoVO&J$CN2jUxv!AA(UtuQe=#gm8*n8>l^FN9r zX}6!ZH(oqqUf$X>;?dC8kJ2xd%(0rWc_XVa8TtI#t-ZU}>M%Z$SL9-OkS=)@G}oB# z{a!oW`jZ}Wn+QGh1Si4hQ&KcYKGhscG6Y5C{AitgqV`mjtK{Zcy1yY}5!%koB-f5F zS@D%IN4-Kh9m!BNPxCp0REmERyF21$7X0atYR*1U1tHGK!JVz*)xtl&PUS8X5MCFd z(pvAozp$apD(ZROMvIotgX~{y;P=GeJlL;erpngfVE<^-94&fsVv5LktS*usY4y^n zKV820m_)I)F0~-C8T&U*Rr@=1i6(3v2#L0?G0NbHiu=&*IgjVXB{A8Z_X%CtQ`9>5 zgV;XxDyjY>ht{P|!y*NiCK7!&Q8hsMl<3 zz7#PT^mc92l6r9aga+p+n?mZt?^#bC3{F?=KYjzt9|`JEAmOt}<7O5kWx7+aiFCAJ z4V?wW1$KlCF&&YHDE*22vFpQKO}RAA;N@V9rQPL@6Zbh$kED|uk?wv0MNrCUq})S; zU(p7$6x%22KQ=Z>=%s_LYs$;@?ec=d*nF&oC#>Gw{jNH-s_>O=G?5a_21$Q%c{j*$Tmy zfFP9;W=V+sKv!1q#cFbCr$BQ5Y2Q4Xq2WVSv)Q&kp{nQe1-O(uv#s^q2FnSN#GOqs zU42L7zqx4wnFJf9`P~9F$*zR?Tz6qn`_1ARN-RdOBF}AnAuF^fh5qck)P<^?a0fTn zS56&V-!t9cK6VCmKS8uSXiTN(eRM#CH_sZ39gM+Z)m^2_zOE5PCbjtK6DHU)FcnA? zR&y60&@wPs>=k^QZTcaI6c|aOx9SM?7(Vy?u6&3jk+%=nGo^(+M9kW>tG*c+tJ>E` z{!xbTq@;A-=^MONR;H-ATYvsiA`;nfbaWHkdax3WcMEVe8-D+es6XCmhm;m&>IFyb zL`NH4t=aR&d-@L=y*l3y-gaxXbHd^&BL-9oqa z_4U;+u6qo-d~R$-(p1_^sS1*(7e`=c7%qY25zSM?t)~e3JdKmFJnyg7^+<{uvFteb zV`G1OfnQNU)7g_Hy3WLzrjvfiMoO4Bvav9RZ2HPw_W6>GmN5wzhho8sRpBpHy{pSg zyb)Y$Yjblc?eCe=RZva*)&Z~EiMJDfK7UQvWLV494FyMNxSM&~8MD|Xv1iina}?>D z5fn|-I?7vfupGo(&}&?@eTigqJs-^Ab)+Y~rA@2p9<5E?%fs71P;VmMlc$nzF`e5_;yHq>sEvDc2;#lQmX6c?N^9zfvcxSa4L!|s{aAq z3Z>H9Gx_PaOP_@LmjfPI4Y&2DAeaz1*1PartLrZ%@BaP|{r?$@`!75TahQ`@2d=ET z!WY_!nZNHpJv{o>K$IER&F6=assNuEd)AjNU4w&-?hc?*gWq(&f0;x5*(U0*P3(?A zi=6|D`Ic-ptgN#Ie1foHiLtNhww!6c3uL22UR3=leS?nJfX)0fGj;?itpmTXPK|CHms?*lFvc67z2$n6Y zy#Vx~|3>rON^Q!Tt02ge_l0*8-SoF`EFRpKn6j#<) zx#|TVXP51tjxHeEXE?YxN^@v;(@eBOU%Xg@aTNu6`)!+bPvA9NQ<);wwwFPY_NjfT z%fUxO-FHqbthZs>Dh#Z4y>E>NuKRP+7mwWM7D}GH@;J{Wo;cpr2HS~^i^n*T*P5|3 zD+8k&p2@6EDeK~;Rz#veTqt+%T3d#ike4|wz3>fDT;yg=x3I~J^^uVa!*aY#75Puc-C&-& zN(!(+e0UzPxlI11X3nQfwNRc8j%T?I?4OT=o`3>SaxW}1_eP0BBFHs_;^ z;B<~DQF;wY#i*?+xOgkhr^mEzw)3zE{xua9-H2FO;xs)6a#j1rVi^_lElwh#O`@jv zXN+LS%>@{nQ~zkA1!w7{JJ*F*#iP~LRS9zVrb{-}>+kf;SwG5Rih7S8Vv>b<@egSC z6zuB5VI=9HLGjAd%ck{gp@gcNCskzK15l}%5$fUfB~ktSn+1_15w zo%(}p;CD>MTzr!P_KD6k;pz&j8+X^+%*sdnem6;s>7}BO{a@py+=zX#*7JP~4t6_) zxs)vh9%G%YlK$HWM`W zUp-pf*)ns|g6+A!Ie(G0{%x7vsDB_XwgOdnhkx|@dZh+}dRd$v(?z>(j|I0tIGHu?_d&@==JHzG^E#a(YQ?BTq3f3o+Tk75idpbes>U6mG76U>EIOD&!hytk8 zy7E)od8cZ$Q7ogBzN(or>z$&V7Bo+tb#gvA{Yw-_tq5k{H;M5JX9oT9u-lA9bYEZD zs-=Tq+A%AwLr?MFCs*-(jrA+ft0}X`EvUFph>8s=zD@Yu8n@S@%WJz^pDm!P^1ALk z4(Nz5Z_Cbqh)Eb!oT``AQ|Mbt-xu-NMG`7_^h0jv`lBcz@2YlCRJvI5qll4`63(`? zx^4?8G;S_Mpr|e6mq1Mim#7ragE88&e6Hy5nXg2s$eT}{*vAVtN#*3e8@!2GT%Sp? zv9Ym@jPhYQA)1_?c3x7l8d0y0OB8e`+pfNA1wiAa_h=%$%+~icXo=#rN6q>ss2pT>0y$YLbDva7bA3K6+?8WmOKyO%E!H#9g=SWFVyoV!ZS7Gnp!O3Yz! zZX;9c#WHXkASDq@ROzU&X3oku(9z*tJHLI{6J-=Gwd2A^`U&}>$l!>JyHCfzx;Ak% z`=Re_}o10#zWq2AZ%(H>Xl8A>#0BO}Dlk!&L)2d4*6tNm6Y zQaPE0&pKuYrxiWtfXM&lY*D1FoA}>YVS81quCx0Xdi5;pwcwRi=g#$sdDM=g-SMH^T8N-Ore>4e zukS5EN5q$@>Vm%3^jS?#y$3f>-npRZh#8+vVGq2IlO(+wi{~F9a*CvcV<{zfBa`{p^3-$atl>M z_>M!duOa*@H_BYQ3={t8(y`~q+Zk?G0mS*$++D>KIh$`Y3`{_|a!alR9nNWxRF)v2 z1e|y%lGJSM!xyKq(dX}1d#K_lzQud#D4t$l@DmHY57tZ&_?I!-tRgMkeuJ) zsN{w%z1B&bG$GdauUNb+XxtTv`uyxt`cdfhSCmwYQ&*E-ot|DF5CeMq2lXju%KYGd z(H1Bk!jrfWq9$QRV~w2&>r!g(a|=;-_~96NCVX4jnxjv*qHy_*H%Y8nfBwPwjJV>T z-@j1^xqT!8N{7gWeH9Jv?md(sNQY9&IkRrxFCnw-;`1wkvD!MaR#&?L=pJlYGo!zb zg%oGi0FR^7aoK8^_u!qpF1uBYFDgpB7zZlk+^revbkh|~NdqFRS;e^jQchJvJFHl5 z*P!^qqt0bxDl&QB;LSvcGF35!SxIe@^IbqhqBr_o5@+gF`rv2?p=vkF>r*+}4uVi9 zeS$kHilq56{~G|2X*!`&8!<+`P3}T^P%pU7TpTP7s4XW)Jp^1vPrJ&8Ow7%&>LA3R z5bXcBkFOqhE!s$AK?bIpg&AL#tJ2`8!&s8)-L_Td)V6x z8*Q#N^L;>WDsm?#Ro!#VA=Tmf;=zs*?$}e{g9GX-QMM%fwdu~A=34@aYU;XbVphiYLx;SiBIbHpd=4za!z0tC zE7F;3glc0xypkBZn&fbw0XIZx>c_@t6+e^R-UL6-Vc(t54=i>i_nE)zPBpbnf7#Ta zlaiM&EPDdF6F$}2)fyZ@h5z)^$3^Tshml5P2{;!tg%?sVHv}QSWB{VLsqb4QeEkNt2 zoZTiKOF}Ffs)6bcDA7Pu(>XG7g>^@15VHqpr86#f_e<(FHDamN68Y`M8?Q-d9a&wD z_WBmevB8F4onb%AP>#a6Y(L}54InoSkX2QTAyQG{FL{~bU(rejX|Q{tNM0g8T-aS6 z;0ChSdj)0#zdbK+Lf@8*DqNjUb&%5GTt&DO4ew`hW~iBpAd`92IdV_fVy^aH+Y&ZE zMXz|?h~v&1P{pj29g?wSUicpP&S)L{E_-424lKTJjllMAHNnHr%k@pXrY{%R2xQ)u zyqD-KK2yyW7eGbgiHUi0y_ZVa0gL(e)o>zM3u3Q0;Ol#i9Q=%c6=(jtp{Qmb0@`1l z&9ixo%g@}VOUi8(*voiA(;8DqfG+y^?zrj&Xy{+C@MULm&iOr7=EZ;O%?fTz`yn;8 z7(En6%gu`>m&w5bR6J5t6o^L7@K@2G=O;1Vrs+Vgxzvr#pX-&2aC|!a$reR6qd2oC zrDB{V6&168!=d!ohH+X=3fR6V1h4(XCeis&P1!i&syR^G!+;dlN0#0{u8C6d| z$*EwUea6C<2d|FWQ{-~7&EF;6ASVj4SE6ZpYjhqQI6fuwT*#9Fl4qcJR~AyJb)&7m zuZXoohgEf9{QoZV*4n5;bjZ*E1BG6vJI@@FIJrDE?4tS_%g}yi5BVax9q~`^r zn~0jK?z_(0K)(*Y*X?4&;|YZ%+c-F^1oPjOu0>l-W}=V^szcQFr<{%~dVVdwvUah3 zVuX>4hC(RdCeZXD;~0}T#GJa+tf<%`AJ~zM5q?9p%o1vD@O8Sh$@Y4uX(Ev#yohL{FVt*(S3p)58R5Wo*$ke#_nvzeEs`rX;P`iJ| z%g@AoCB{V~{(m0F)-3-t`AQ0@nda-s{;v98Ej~|% ztHzMx@wnFpNy<^ifYVk^uQ2&3w%Y^m9E#ltB>fz#-EO=BizQO>&7PE zS3`bxUBc6TW`HfoF-c}9l{sx@F?ip`QEe0buxB+a-4i6_a%6k^T=dD~=wFjnrcth{ zvLgon$^y^k4 zvL)m-MlHMedQKBo!{!*JloY+*!}>n?<%Kr@mmV~{^4M=47|gfkzW*C#k(@}uUf2%^ zwsSxB61UFkRUFz_Kw~iJuQ6l;CV}_;kCT6H`OWAns@!aWn}QHAY?GT>037puaR13; z@<4kfg1H-SBEkGt1mXYo^3;mqkuKaJrsp&WFEgj&M_WlzE`f#xq1?ImI&sXtUz`C;UM4iQiymEQ^>A$!zqoOpsRKNJ_QqLL}Dqx%uo(h}2_2E1Q(YQ?gX zZRQVQd%ww`S4(4A{Bt+g1gsu(c4W5JKG3N>SKD4d5G@Eb4?PMCt$V)VJz!#@%Pgu4 z93EByCdk)>ZQ4eT2i;UizJIRWJ9Y1v25ERQ?l!7X{bt>jJaZk zoY%$BWA{zdb&GJ>Tc0`8bv%4rqDO1o@r4a2Zv!5&TlBv7JUuF|Y44hVVmu}{Io5Ud zruug1tgf$*aMVEi*FPyKo?6OfGxOJ8&h_j?;#>kB)#3r}QBjMD-g{av0Av6?uTzta zz^*`a(hoY*ZYM%{IWT!uSo#h+W1!~(b6YlNM>v`1S z-=m?U`K}9m@Ae<`j5$Bbo&oUMl5;0$u*oa|knmFyW>*w`9^MSK0;Or|vMN=GH6$GN zHb5k!t)&;v253ICyGB;tA;q z!TOl($O*=~Do!960o1JT%KK0Xo828H>mB8MwF!Q{dkdqvnWnvSE%%6tae+l#z17X* z^yb5T$%KM(4Hezq_5CKToSi7WdREg@`G^z&T=#)f^pfP9j!53_e8$u7Cm1Nt#|z6n zz7!VUS9INNc4x?KFZtrP@#YI8OIuw9#es$OLl^rxfp33b7f3veTcF0|!#n5Z(y)V7 zBFP&3F2$^RHf1ZtN#@djgMAVGJ&L+Il|-B$J+_0E2QjwWaW@eR?~G7R|k4Q^>ciTA7`lj#mrWoOi) zlX?bR9UR0gb7$p1FF^~{8~om)>vfbGh)(Vy4tcC&HIrDFVlR9VGN7S(Z1`@(;Tc75 zUnW)qXbKo@9gc^8C|pcUf7sa2&@--1yOKLDIncppJ@KL5H^~j(%o0GH1a!V-Y#f{S zy(2y0z;1Wp(~B$TMz535wk^*XE8v^|{tQWAs-^>4Obs#)6BARJDAS`BcCo+5)%|9M%L_kJ~&nRdwP0$A1|1cG?!{tKQx7QO4;!jutb zTz59r)<(PW(2KOmZZ^r8h-9t@AQsiLR~JzBJnP@ZhmpkdZQv1uUf=K&l~#=2@VPbW z^=T`W8RnA=)r7r34C1$jZH5|5T4Eb9i!K|(pLoxNJM~cURb!H6EiI*;HQ_rdu}e?L zS?jy;lR-=fjf_35wlnojYw2hGQffOzN>#hYs=B#&bso#QkeHSd`@2K z9g2{jtd@kFDk`@9R7iWDb0d+LcqjS44Y;?=s%k)(3)A;mmAZa1z+AEGS#pP#0lbF}l!$j|BWa5~&sY zf0dCA!Ex8-dk9s&>D>`tfJ7AXch#D`#+yF8YM-;nj3StC37xrUIE)^|yOOKO)jUAQ zdX718%EQ;Nu(1;uP{B?VDcDM_FZw6O$Tg~&U2OJ$K=T-XTos#>M#?r}*0RtO+Sl;Zp1#X}CuP)vu}U3X&n zEx9^)oSV&LrzUxmtlQcJ-~L(36}d^@nnB3)2_}BO+4qsnbC2amcHe4#IVgw%gbKJV#%c>qvOcXA?M--Jjy9nmP9K z^1>7c7Z1P*Qj*E^-~9%DCC}NF{ALodymG%fdt6KJJTVf&8#2~R0Ct^*z`cSt}XGRkh}*^{hm+t%~zFZY$EqJLFhF!Fkl};-9MeFM&pZu@}BRl%QYn zP5rtw7aZoaqx}mBvAD4dPQLV$AlF`|JT*F?r9XAv^DCP}{w@!N7dBr#fi`8fXO?I* zgEfipR_on>twt(ojt_tg2q%qyKhwH#9`iHgnG{*9#W0fMn59HhdwUwzg+*vdP21Vi+BKQ;mi56}sNfGm zKaWXXa5Pb4#mDSrwXQn_I&OEVZ9CUJe~#$Wy@3Zaftz40ze2=02XYt~lx?wqPdnM5 zA+AF_oE{GQ>0U<7YE`<(3%$8=zQa~jC#3P^EcUqIRar9LSes35flxnVwXzB5iAk(u zbG0MeABq>cy_rmZrnr~GSoD>P4xm@m8iCrfN%MR5moAK;9<`2#t%f zR#sQQ_8nQ-t&co*MlxID?jG0m41mW2bPU#=vLVpX0lH4WzX=8w(i`CJtTB3gE182P zALu0pgEg==VueS8_dA8jUMqAOR~nGuAV`9}$UZTSNb|i95)BQWR~M|_OE_zcdJ4qc zXN6imSq0I0SZl5q8tSy=%Ky+~bDp^jXU(99;tA(Uzlk>j{2qji=D&WuJtik;$tJd| z{8nCqiKInnq|1-keVJu@@) z?uJq8!=&sVQ1x*MXZ!Uo5EMqAJrMM{Ho@|L_PYS&s2L=Dz}x?RL}XTU85}6n6ZXBa zC?X5xjM%J*YG+Oy1ejY7z}2!{PHmnee_ABQu^G-i5qUiH1qcb>-hCtxJVQvNF8zxrk5) zaM|bj!5>tNI*v(-+Hu=ckO{B)($@bQ(LDF&tJz;Lr0n@z-gmfcW91<75}0yl$&YB; z5Ijc?%B)o8G*vIO7}Pd=9cBc|ymjl059B^q>%+PzZte$^T$;rE=Y)W4 zZ1dVGsuXnx%;B(;qLmed{C)>TvM{`ut?@(Rs;-vTrv@?mE(EVeOD9`V|E z^xdOKqHuc$8X~VNdM5Y;)Yoy*`z0*mwf8wew<{m}Xp@nQOy!DxuTh`u0NYW?Rg&Yn zZbiJ`#r@$fM`g;r*;Zf&w;sTlUN8neTQUF&N{p*L6c)H`&T4$9hKTV2Bz>`ys$O(` zY@MZ}HwOF9mOq1QVWDYz%TRu|w&D2bPY~tIdvW$s*!O}E@fTe3G4ZoD{|}4DhbMWG z^S+VV=e3v65*!lN2=Kg=7p2u*NTO-^?hPMPJTL--mJt{PIP19Y>R-P@txa1tg}9-0bOW)quEc zG`eCL7)=HlPUL}{@xAbl8FEJW+TDLOeT+R>mU{onZdeeRkI&;Uf(W1n`%+|p*pu$i zHW>N%kuGAu%h!lO*zW1nvn5y@)^$4v(%RWAxSkf{96$ThAqQdV$C&8$e~xqtaHECg z@-=$Cz%?KI`HvPnI{SaJ;5nmiDi={tsoNY){GNifbK(Bv$)_*Xq$#Kg!NpeY*Ix(d`;!BBQAxqijU6aiV@|XG=#%OuR!%0&*QZ`Z8oB zDJeW0TNo1%2N(78RJbr#oV5j`0D!5kxeD#$^ zyW7mpJg>X{cFboJn5)NuKYbWD^547$VGcEN$RQhmBl<}|AqM8=Cl{2+g>X`lfYJrnhYVf} zp2%ZVh;pc&mL)9^8F)O-pVj|lQ{&>=aNz@HeCOhd7{cW#oWY_IK6`1dK+-`?Gw~-e zF99Ez_!}_G07CXSBJZ8IzYA%+YO#%(O@M2KXFEq`wrhvDh) zb6RSxjrJ}p-nj4_h?#@WGeHwh@sw(*XyCST&$H!5ms&kkXkXwELfU-ZBHxai2v|oWFuTd`eUkL(Q z{fm{JwUUy_9IDqw5Nlaxa=%u{f5#fvM0*Sx>Gu1G0OxD9F(qt(!$w`3JLHo;fIZXn zZ00SnUy;6ARC(X_6N^QVg)6>sfSPigI&*?eSkvtuh}Ln?hHvm+>;}vCl3|JDG0i?9 z<@)$5q@_srXUHH$DX9()9X#)m-fEhT0N%CwK;SG?25Y#QHt%U%~VS zREC3wtrtmzn;=sq>D_o}8ryAwo=y)ovdI5P=q>{OA#}ueRX0nsHRr~CI%q_^)x3q1 zCP1!W?dnV{tP4+D0O49TyUU1223f+s*J(|+4<1Lb2ra85=E`E)1JDog^OlQ*Sb>=u zb9|E26hy2Y(OdsfhN|KAkCUR5*dYuboJ8n|Iz7=LzM@Q z)4~Isp+bg|n$}47-Z0XI44?}!oXLI^c5q=`p5!4Cb&dugh9S~*;LEpe8dMpr5~DnM zKxEz1yQT#2Jp!wEe8u`0K29ETx$Ta};8#8t5nxS+U%|BCY57vh7143S!U;jn9tZqH`fN^K3{=7`MSBB$5o0k5u9!9W6fszS#AF*m~c z4O-_ocE$x>d3y7aL>wpTr)MIS41WpKe? zA|N10%S2b?s#CHB?*s}9na?I?QqII1;3igP1xrVSzI;s;<}>7FeG?rRkcADO_2nM1 zZ>X3%*#=hb|HYfCimIV5xKEMc^ZKk>FJJ9aB71o7@bG1*i)76qMPi*<2;FIpF5zs}VNTu?fgVE;~7Q&#i!DyU2_mWA7`#|(# zq@olPQ|N}|=Fvk=xDE?(C?Ti?{bI8ZOIroe*jSuhoNOg9s^fCjpj4KBVaaQ;k1vCU z0XYN_j>!0to{og%b>X(iGoitia0&xPEOsyznOd}5^YxjUoC?oR;iRv%1WWGwMdeR? zz`8p}RNY2?G3 z6B(A~wcV|T6w zD!;6ZE`wtNtRn|JB&ihy>vcGn`gd=z5SW9Y_aug2o8k0P{kx;fbx01xB6bpI>**+p zW$GmmtOcht*L3?cJO~_L8B=-r8IEm0Tb&t&nAun4<88i(g0q7^GURW zUUER>K^29>O9PTb39?l{T4i-9rb0M=3e}76gZuReeRhwD!tWvC#hC{HHnV`>K^d46 zLVJvw=at@slp^Ga76~H*Lnm@{uUI`XeLrN-*s95EjnM(qhJX&M4va4@5CdEFPl(xk zlq|TzWfBOCtJ)EL2;x_dN%=lT1*js@9I?Q+UIXyk&>)nax^=EtG zA{h(4q025iQ+-?JqImz)z-eT`NQ=)EnvU~W6+xN35{p&jS~t^2#B(o6Y^QWv2e~QB zR($W9A+N^DQ|H^6%U2-w-d}Jszc|-l0fx}B-EJ@0b%YFpVRIjf4eSkpJ@0t(=LIcQ zScks8DGwMzlgd$9#3pEqKS^ix+dTB#j`@e1g4)Ke;s26#6-&#^i~@{6#6!zb+Mw6; zFBX+3Pn-GasaR1TsF!r-|#PTmR-6IbNzIza#8dZaf=A)BzsvM9WZU(}PF4R*|vo{-BN?J6qeN z{gdwZQp|0?sXoOdm2$FCiws-^vc0sjG9rx0K2~H2B9n%Py}Uu`&dxsp<}3a@`eF0m zOoOI5cFft1jBGKa8G$JxH&=jf;v7AouRTRQ)M{R{1?U3_ zM7@AebY1W9>R+O2a$USUP#vBS>*#o1oD-8BogRp0D8+!~Hdt%@H7PBsCc1P8P6D*C z<N-f^dc z`LWLaZc3Ml5TbQ3D*GE6o^G-l8?Br<_RTwTSTXIzR_VM}*|oUz7GeQwv}IB$s+n7x zcO2{Nb6k6aOwg~#+XU8m;}#Al7CklFLS_f__M`PFA%0Q>WNB$Qz{AvgI5EL_h>5;j zKfNEHR19v07(E1$wn03^;!t(>c0gFcvFWM*2cs4L4}29gX6|5<=Ua%xV_zyJULGXmL6K+FfkHXuRW)ST4Z s)VvZqpa4)UCy*#Ej4v)J%FIiLX#g3JmQ#{Vk|uVbru4khf}H#k0D5999{>OV literal 0 HcmV?d00001