diff --git a/lecture_15/lecture_15.pdf b/lecture_15/lecture_15.pdf
new file mode 100644
index 0000000..94d5f8c
Binary files /dev/null and b/lecture_15/lecture_15.pdf differ
diff --git a/lecture_16/.compiled_data.csv.swp b/lecture_16/.compiled_data.csv.swp
new file mode 100644
index 0000000..70ab7de
Binary files /dev/null and b/lecture_16/.compiled_data.csv.swp differ
diff --git a/lecture_16/MonteCarloPi.gif b/lecture_16/MonteCarloPi.gif
new file mode 100644
index 0000000..8ffb745
Binary files /dev/null and b/lecture_16/MonteCarloPi.gif differ
diff --git a/lecture_16/MonteCarloPi_rand.gif b/lecture_16/MonteCarloPi_rand.gif
new file mode 100644
index 0000000..aaa31d7
Binary files /dev/null and b/lecture_16/MonteCarloPi_rand.gif differ
diff --git a/lecture_16/gen_examples.m b/lecture_16/gen_examples.m
new file mode 100644
index 0000000..c5d6890
--- /dev/null
+++ b/lecture_16/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_16/lecture_16.ipynb b/lecture_16/lecture_16.ipynb
index 74bdda2..be266b3 100644
--- a/lecture_16/lecture_16.ipynb
+++ b/lecture_16/lecture_16.ipynb
@@ -1,10 +1,86 @@
{
"cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Curve-Fitting\n",
+ "\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": [
+ "# Statistics and Curve-Fitting\n",
"## Linear Regression\n",
"\n",
"Often, we have a model with empirical parameters. (e.g. Young's modulus, Poisson's ratio, drag coefficient, coefficient of restitution, spring constant)\n",
@@ -17,7 +93,9089 @@
" \n",
"- Factors not accounted in model (e.g. 2D effects of 1D approximation)\n",
"\n",
- "These can lead to **noise** (lack of precision) and **bias** (lack of accuracy)"
+ "These can lead to **noise** (lack of precision) and **bias** (lack of accuracy)\n",
+ "\n",
+ "Consider a piece of glass being stretched. \n",
+ "\n",
+ "\n",
+ "\n",
+ "It is clear that a straight line is a \"good\" fit, but how good and what line?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Statistics\n",
+ "\n",
+ "How do we describe *precision* and *accuracy*?\n",
+ "\n",
+ "- mean\n",
+ "\n",
+ "- standard deviation\n",
+ "\n",
+ "Take our class dart problem\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "darts=dlmread('compiled_data.csv',',');\n",
+ "x_darts=darts(:,1).*cosd(darts(:,2));\n",
+ "y_darts=darts(:,1).*sind(darts(:,2));\n",
+ "\n",
+ "colormap(colorcube(length(darts(:,3))))\n",
+ "\n",
+ "scatter(x_darts, y_darts, [], darts(:,3))\n",
+ "xlabel('x position (cm)')\n",
+ "ylabel('y position (cm)')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "mu_x = 0.90447\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "mu_x=mean(x_darts)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "s_x = 4.2747\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "s_x=std(x_darts)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "mu_y = 0.88450\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "mu_y=mean(y_darts)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "s_y = 4.6834\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "s_y=std(y_darts)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 83,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x_vals=linspace(-15,20,30);\n",
+ "hist(x_darts,x_vals,1);\n",
+ "[histFreq, histXout] = hist(x_darts, 30);\n",
+ "binWidth = histXout(2)-histXout(1);\n",
+ "bar(histXout, histFreq/binWidth/sum(histFreq));\n",
+ "pdfnorm = @(x) 1/sqrt(2*s_x^2*pi).*exp(-(x-mu_x).^2/2/s_x^2);\n",
+ "%cdfnorm = @(x) 1/2*(1+erf((x-mu_x)./sqrt(2*s_x^2)));\n",
+ "%hist(x_darts,x_vals,trapz(x,f))%,cdfnorm(max(x_darts))/2)\n",
+ "hold on;\n",
+ "plot(x_vals,pdfnorm(x_vals))\n",
+ "n=2; % n=1, 68% confidence, n=2, 95% confidence, n=3, 99% conf\n",
+ " plot([mu_x+n*s_x,mu_x+n*s_x],[0,0.1],'r-')\n",
+ " plot([mu_x-n*s_x,mu_x-n*s_x],[0,0.1],'r-')\n",
+ "\n",
+ "xlabel('x position (cm)')\n",
+ "ylabel('relative counts')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 84,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "y_vals=linspace(-15,20,30);\n",
+ "hist(y_darts,y_vals,1);\n",
+ "[histFreq, histXout] = hist(y_darts, 30);\n",
+ "binWidth = histXout(2)-histXout(1);\n",
+ "bar(histXout, histFreq/binWidth/sum(histFreq));\n",
+ "pdfnorm = @(x) 1/sqrt(2*s_y^2*pi).*exp(-(x-mu_y).^2/2/s_y^2);\n",
+ "%cdfnorm = @(x) 1/2*(1+erf((x-mu_x)./sqrt(2*s_x^2)));\n",
+ "%hist(x_darts,x_vals,trapz(x,f))%,cdfnorm(max(x_darts))/2)\n",
+ "hold on;\n",
+ "plot(y_vals,pdfnorm(y_vals))\n",
+ "n=2; % n=1, 68% confidence, n=2, 95% confidence, n=3, 99% conf\n",
+ " plot([mu_y+n*s_y,mu_y+n*s_y],[0,0.1],'r-')\n",
+ " plot([mu_y-n*s_y,mu_y-n*s_y],[0,0.1],'r-')\n",
+ "\n",
+ "xlabel('x position (cm)')\n",
+ "ylabel('relative counts')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 76,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x_exp=empirical_cdf(x_vals,x_darts);\n",
+ "plot(x_vals,x_exp)\n",
+ "hold on;\n",
+ "plot(x_vals,normcdf(x_vals,mu_x,s_x),'k-')\n",
+ "legend('experimental CDF','Normal CDF','Location','SouthEast')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 170,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% plot the distribution in x- and y-directions\n",
+ "gauss2d = @(x,y) exp(-((x-mu_x).^2/2/s_x^2+(y-mu_y).^2/2/s_y^2));\n",
+ "\n",
+ "x=linspace(-20,20,31);\n",
+ "y=linspace(-20,20,31);\n",
+ "scatter3(x_darts, y_darts,ones(length(x_darts),1))\n",
+ "xlabel('x position (cm)')\n",
+ "ylabel('y position (cm)')\n",
+ "hold on\n",
+ "[X,Y]=meshgrid(x,y);\n",
+ "view([1,1,1])\n",
+ "\n",
+ "surf(X,Y,gauss2d(X,Y))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Monte Carlo Simulations\n",
+ "\n",
+ "Monte Carlo models use random numbers to either understand statistics or generate a solution. \n",
+ "\n",
+ "### Example 1:\n",
+ "#### Calculate $\\pi$ with random numbers. \n",
+ "\n",
+ "Assuming we can actually generate random numbers (a topic of philosophical and heated debates) we can populate a unit square with random points and determine the ratio of points inside and outside of a circle.\n",
+ "\n",
+ "\n",
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The ratio of the area of the circle to the square is:\n",
+ "\n",
+ "$\\frac{\\pi r^{2}}{4r^{2}}=\\frac{\\pi}{4}$\n",
+ "\n",
+ "So if we know the fraction of random points that are within the unit circle, then we can calculate $\\pi$\n",
+ "\n",
+ "(number of points in circle)/(total number of points)=$\\pi/4$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 96,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "function our_pi=montecarlopi(N)\n",
+ " % Create random x-y-coordinates\n",
+ "\n",
+ " x=rand(N,1);\n",
+ " y=rand(N,1);\n",
+ " R=sqrt(x.^2+y.^2); % compute radius\n",
+ " num_in_circle=sum(R<1);\n",
+ " total_num_pts =length(R);\n",
+ " our_pi = 4*num_in_circle/total_num_pts;\n",
+ "end\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 98,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "mean value for pi = 3.132400\n",
+ "standard deviation is 0.054925\n"
+ ]
+ }
+ ],
+ "source": [
+ "test_pi=zeros(10,1);\n",
+ "for i=1:10\n",
+ " test_pi(i)=montecarlopi(1000);\n",
+ "end\n",
+ "fprintf('mean value for pi = %f\\n',mean(test_pi))\n",
+ "fprintf('standard deviation is %f',std(test_pi))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example 2\n",
+ "#### Determine uncertainty in failure stress based on geometry\n",
+ "\n",
+ "In this example, we know that a steel bar will break under 940 MPa tensile stress. The bar is 1 mm by 2 mm with a tolerance of 1 %. What is the range of tensile loads that can be safely applied to the beam?\n",
+ "\n",
+ "$\\sigma_{UTS}=\\frac{F_{fail}}{wh}$\n",
+ "\n",
+ "$F_{fail}=\\sigma_{UTS}wh$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 112,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=1000;\n",
+ "r=rand(N,1);\n",
+ "wmean=1; % in mm\n",
+ "wmin=wmean-wmean*0.01;\n",
+ "wmax=wmean+wmean*0.01;\n",
+ "hmean=2; % in mm\n",
+ "hmin=hmean-hmean*0.01;\n",
+ "hmax=hmean+hmean*0.01;\n",
+ "\n",
+ "wrand=wmin+(wmax-wmin)*r;\n",
+ "hrand=hmin+(hmax-hmin)*r;\n",
+ "\n",
+ "uts=940; % in N/mm^2=MPa\n",
+ "\n",
+ "Ffail=uts.*wrand.*hrand*1e-3; % force in kN\n",
+ "hist(Ffail,20,1)\n",
+ "xlabel('failure load (kN)')\n",
+ "ylabel('relative counts')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Normally, the tolerance is not a maximum/minimum specification, but instead a normal distribution that describes the standard deviation, or the 68 % confidence interval.\n",
+ "\n",
+ "So instead, we should generate normally distributed dimensions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 113,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=1000;\n",
+ "wmean=1; % in mm\n",
+ "wstd=wmean*0.01; % standard deviation in mm\n",
+ "hmean=2; % in mm\n",
+ "hstd=hmean*0.01; % standard deviation in mm\n",
+ "\n",
+ "\n",
+ "wrand=normrnd(wmean,wstd,[N,1]);\n",
+ "hrand=normrnd(hmean,hstd,[N,1]);\n",
+ "uts=940; % in N/mm^2=MPa\n",
+ "\n",
+ "Ffail=uts.*wrand.*hrand*1e-3; % force in kN\n",
+ "hist(Ffail,20,1)\n",
+ "xlabel('failure load (kN)')\n",
+ "ylabel('relative counts')"
+ ]
+ },
+ {
+ "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": 123,
+ "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",
+ "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)"
]
},
{
diff --git a/lecture_16/q1.png b/lecture_16/q1.png
new file mode 100644
index 0000000..f086fc8
Binary files /dev/null and b/lecture_16/q1.png differ
diff --git a/lecture_16/q2.png b/lecture_16/q2.png
new file mode 100644
index 0000000..6f91d19
Binary files /dev/null and b/lecture_16/q2.png differ
diff --git a/lecture_16/sgs_strain.avi b/lecture_16/sgs_strain.avi
new file mode 100755
index 0000000..625d925
Binary files /dev/null and b/lecture_16/sgs_strain.avi differ
diff --git a/lecture_16/sgs_strain.gif b/lecture_16/sgs_strain.gif
new file mode 100644
index 0000000..8fa0038
Binary files /dev/null and b/lecture_16/sgs_strain.gif differ