{"id":141,"date":"2008-06-11T10:56:52","date_gmt":"2008-06-11T15:56:52","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/11\/interpolation-in-matlab\/"},"modified":"2018-01-08T14:46:18","modified_gmt":"2018-01-08T19:46:18","slug":"interpolation-in-matlab","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/11\/interpolation-in-matlab\/","title":{"rendered":"Interpolation in MATLAB"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>I'd like to introduce a new guest blogger - <tt>John D'Errico<\/tt> - an applied mathematician, now retired from Eastman Kodak, where he used MATLAB for over 20 years. Since then, MATLAB is\r\n         still in his blood, so you will often find him answering questions on the <a><tt>newsgroup<\/tt><\/a> and writing new utilities to add to <a title=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadCategory.do (link no longer works)\">MATLAB Central<\/a>.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#2\">Polynomial regression<\/a><\/li>\r\n         <li><a href=\"#3\">Plot your data<\/a><\/li>\r\n         <li><a href=\"#6\">Look at the residuals<\/a><\/li>\r\n         <li><a href=\"#9\">Higher order polynomials - are more terms always better?<\/a><\/li>\r\n         <li><a href=\"#13\">What do you do with interpolation?<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>I'll assume you have some data points through which you wish to pass a curve, interpolating your data. (Initially, I will\r\n      only talk about problems with one independent variable.) In these coming blogs, I'll try to show some ways to do exactly this,\r\n      i.e., find a curve that passes through your data. Along the way I'll try to give some pointers on curve fitting, interpolation,\r\n      modeling, approximation, etc.\r\n   <\/p>\r\n   <h3>Polynomial regression<a name=\"2\"><\/a><\/h3>\r\n   <p>A valid question for some to ask is why start out with a discussion about polynomial <a href=\"http:\/\/en.wikipedia.org\/wiki\/Regression_analysis\"><tt>regression<\/tt><\/a> , when we really wanted to talk about interpolation. Many people mistake the ideas of interpolation with the approximation\r\n      produced by a regression model, calling both of these things interpolation. So I'm starting out with some discussion about\r\n      what interpolation is not. Plus, I want to assure an understanding of polynomials, since many of the tools for interpolation\r\n      are polynomial based in some way.\r\n   <\/p>\r\n   <p>Let us start by creating some data. An exponential is a good place to start, a simple curve shape that is easy to fit.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">x = -1:.1:1;\r\ny = exp(x);<\/pre><h3>Plot your data<a name=\"3\"><\/a><\/h3>\r\n   <p>It is always a good idea to plot your data. In fact, I'll suggest that you should plot everything. Plots are useful, since\r\n      your eye and brain are splendid at things like pattern recognition. Only you know your data, as the scientist, engineer, or\r\n      analyst. You will always benefit if you can employ your knowledge of a system as part of the modeling process.\r\n   <\/p>\r\n   <p>Next, always think in advance about your goals for any model.<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>Will you use this model purely for interpolation, i.e., for predictive purposes only?<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>Do you need to derive some understanding about the process from your model? Perhaps you need to estimate an upper asymptote\r\n            of your process. If so, then you may want a model that has such an upper asymptote built into it.\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>Will you need to include the model coefficients in some paper that you expect to write? Will you need to use this model in\r\n            MATLAB, or in some other tool?\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>Must the model be simple to evaluate?<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>Must the model be efficient to evaluate? You should know whether you will need to evaluate this model millions of times, or\r\n            just once.\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>What do you know about your data? Is there noise in the data? May you ignore that noise, or must you smooth the noise away?<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>What do you know about the underlying functional relationship? Is it monotone? Increasing? Decreasing? Positively\/negatively\r\n            curved? Must it pass through a specific point?\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>Must the interpolant have specific requirements in terms of continuity or differentiability?<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>For example, I once had a problem where I knew I had some significant noise in our process, but I chose not to do any smoothing\r\n      anyway. Any such smoothing would also have smoothed out some potentially important features of our process. Since I could\r\n      survive with the noise in my interpolant, I chose the lesser evil.\r\n   <\/p>\r\n   <p>Always know your goals for any such task.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">plot(x,y,<span style=\"color: #A020F0\">'bo'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Y<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Exponential data'<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_01.png\"> <p>This is a nice, well-behaved function. It is of the form y = f(x). I'll pretend for the moment that I have no idea what was\r\n      in the underlying functional relationship.\r\n   <\/p>\r\n   <p>One thing I learned in some long past calculus course was that a Taylor series will provide an approximation for many functions.\r\n      Polynomials are useful things. They are simple to use, simple to build, simple to work with. And a truncated Taylor series\r\n      is basically a nice, simple polynomial.\r\n   <\/p>\r\n   <p>So we will start with a linear polynomial approximation for this curve, built using <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/polyfit.html\"><tt>polyfit<\/tt><\/a> . This is a utility provided in MATLAB to estimate a polynomial model using linear regression techniques. We could also use\r\n      many other tools to build our polynomial model, but polyfit is a useful one, and easy to use.\r\n   <\/p>\r\n   <p>A linear, or first <a href=\"http:\/\/en.wikipedia.org\/wiki\/Degree_of_a_polynomial\">degree<\/a> polynomial (many use the words \"order\" and \"degree\" interchangeably), might be written mathematically as <tt>y(x) = a1*x + a2<\/tt>. In MATLAB we will merely store the coefficients, as a vector [a1,a0]. Note that a polynomial in MATLAB has it's coefficients\r\n      stored with the highest order term first.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">P1 = polyfit(x,y,1)<\/pre><pre style=\"font-style:oblique\">P1 =\r\n    1.1140    1.1937\r\n<\/pre><p>We can evaluate the polynomial with <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/polyval.html\"><tt>polyval<\/tt><\/a>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">yhat = polyval(P1,x);\r\n\r\nplot(x,y,<span style=\"color: #A020F0\">'bo'<\/span>,x,yhat,<span style=\"color: #A020F0\">'r-'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Y<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Linear polynomial fit'<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_02.png\"> <h3>Look at the residuals<a name=\"6\"><\/a><\/h3>\r\n   <p>In fact, I'll claim the relationship we are modeling is not terribly well represented by a linear model. Depending on your\r\n      needs for this model, you might have decided differently.\r\n   <\/p>\r\n   <p>When you build a regression model, look at the residuals.<\/p><pre>ALWAYS PLOT EVERYTHING!<\/pre><p>Plot your residuals. In general, some good ways to plot the residuals are versus<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>the independent variable. Look for patterns. Patterns here in the residuals are often a hint that you should look more deeply.<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>the measurement order, just in case there was a problem with your equipment. (I've seen many cases where measurement apparatus\r\n            was re-calibrated at the end of every day, but an experiment spanned more than one day.)\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <div>\r\n      <ul>\r\n         <li>the dependent variable. This might help pick out cases of non-uniform variance.<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>Look at your plots. Think about what you see there. Compare it to your expectations.<\/p>\r\n   <p>Since this data was very simply generated, I'll dispense with some of those plots for brevity. Note that the residuals for\r\n      this linear fit look vaguely like a quadratic polynomial. This is often the case when there is lack of fit in a polynomial.\r\n      That lack of fit often looks like the first term we truncated from the Taylor series.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">res = y - yhat;\r\nplot(x,res,<span style=\"color: #A020F0\">'bo'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Residuals<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Residuals for the linear fit'<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_03.png\"> <p>If the residuals looked vaguely parabolic in shape, then it might make sense to use a second order (quadratic) polynomial\r\n      for our fit.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">P2 = polyfit(x,y,2)\r\nyhat = polyval(P2,x);\r\nplot(x,y,<span style=\"color: #A020F0\">'bo'<\/span>,x,yhat,<span style=\"color: #A020F0\">'r-'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Y<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Quadratic polynomial fit'<\/span><\/pre><pre style=\"font-style:oblique\">P2 =\r\n    0.5402    1.1140    0.9956\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_04.png\"> <p>Note that the residuals here look vaguely like a cubic polynomial, although they are much smaller in magnitude than the previous\r\n      fit. Again, at each step as we increase the order of the model, the residuals will often tend to look much like a polynomial\r\n      of the next higher order.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">res = y - yhat;\r\nplot(x,res,<span style=\"color: #A020F0\">'bo'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Residuals<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Residuals for the quadratic fit'<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_05.png\"> <h3>Higher order polynomials - are more terms always better?<a name=\"9\"><\/a><\/h3>\r\n   <p>Polynomial modeling with polyfit is indeed simple and easy to do. In fact, we might decide to just use higher and higher order\r\n      polynomials, always chasing the term we truncated in the previous model. After all, if a quadratic model is better than a\r\n      linear one, then why not go to a cubic? Ten terms must be better than nine. Look at a tenth order model.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">P10 = polyfit(x,y,10);\r\ndisp(P10(1:6))\r\ndisp(P10(7:11))<\/pre><pre style=\"font-style:oblique\">    0.0000    0.0000    0.0000    0.0002    0.0014    0.0083\r\n    0.0417    0.1667    0.5000    1.0000    1.0000\r\n<\/pre><p>We can write the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Taylor_series\"><tt>Maclaurin series<\/tt><\/a> representation for the exponential function as\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_eq98064.png\"> <\/p>\r\n   <p>We can compare P10 to the coefficients of the known series. How well did we do in the fit?<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">series = 1.\/factorial(10:-1:0);\r\ndisp(series(1:6))\r\ndisp(series(7:11))<\/pre><pre style=\"font-style:oblique\">    0.0000    0.0000    0.0000    0.0002    0.0014    0.0083\r\n    0.0417    0.1667    0.5000    1.0000    1.0000\r\n<\/pre><p>Note that the higher order coefficients deviate somewhat from the known series, although the lower order terms appear to be\r\n      quite accurate.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">yhat = polyval(P10,x);\r\nplot(x,y,<span style=\"color: #A020F0\">'bo'<\/span>,x,yhat,<span style=\"color: #A020F0\">'r-'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Y<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Tenth order polynomial fit'<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_06.png\"> <p>The residuals oscillate tightly around zero. In fact, they are so small that this last polynomial begins to approach a true\r\n      interpolating polynomial. Perhaps if we just add a few more terms we may get there. The numerical issues of floating point\r\n      arithmetic will often preclude true interpolation down to the least significant bit anyway.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">res = y - yhat;\r\nplot(x,res,<span style=\"color: #A020F0\">'bo'<\/span>)\r\nxlabel <span style=\"color: #A020F0\">X<\/span>\r\nylabel <span style=\"color: #A020F0\">Residuals<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\ntitle <span style=\"color: #A020F0\">'Residuals for the tenth order fit'<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/141\/blogJD1_07.png\"> <h3>What do you do with interpolation?<a name=\"13\"><\/a><\/h3>\r\n   <p>I'll start talking about true interpolation in my next blog. But remember that interpolation is different from the approximations\r\n      provided by polyfit or any other regression modeling tool.\r\n   <\/p>\r\n   <p>Until that time, please give me your <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=141#respond\"><tt>comments<\/tt><\/a> on this blog, or ideas for future blog topics on interpolation or modeling in general. Do you have some interesting applications\r\n      of interpolation? Some interesting problems?\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_90e4d2c2328643a19d71a40dc659db85() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='90e4d2c2328643a19d71a40dc659db85 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 90e4d2c2328643a19d71a40dc659db85';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Loren Shure';\r\n        copyright = 'Copyright 2008 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_90e4d2c2328643a19d71a40dc659db85()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.6<br><\/p>\r\n<\/div>\r\n<!--\r\n90e4d2c2328643a19d71a40dc659db85 ##### SOURCE BEGIN #####\r\n%% Interpolation in MATLAB\r\n% I'd like to introduce a new guest blogger -\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadAuthor.do?objectId=801347&objectType=author |John D'Errico|>\r\n% - an applied mathematician, now retired from Eastman Kodak, where\r\n% he used MATLAB for over 20 years. Since then, MATLAB is still in his\r\n% blood, so you will often find him answering questions on the\r\n% <http:\/\/ |newsgroup|>\r\n% and writing new utilities to add to <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadCategory.do MATLAB Central>.\r\n\r\n%%\r\n% I'll assume you have some data points through which you wish to\r\n% pass a curve, interpolating your data. (Initially, I will only\r\n% talk about problems with one independent variable.) In these\r\n% coming blogs, I'll try to show some ways to do exactly this, i.e.,\r\n% find a curve that passes through your data. Along the way I'll\r\n% try to give some pointers on curve fitting, interpolation, modeling,\r\n% approximation, etc.\r\n\r\n%% Polynomial regression\r\n% A valid question for some to ask is why start out with\r\n% a discussion about polynomial\r\n% <http:\/\/en.wikipedia.org\/wiki\/Regression_analysis |regression|>\r\n% , when we really wanted\r\n% to talk about interpolation. Many people mistake the ideas of\r\n% interpolation with the approximation produced by a regression\r\n% model, calling both of these things interpolation. So I'm starting\r\n% out with some discussion about what interpolation is not. Plus, I\r\n% want to assure an understanding of polynomials, since many of the tools\r\n% for interpolation are polynomial based in some way.\r\n%\r\n% Let us start by creating some data. An exponential is a good place\r\n% to start, a simple curve shape that is easy to fit.\r\n\r\nx = -1:.1:1;\r\ny = exp(x);\r\n\r\n%% Plot your data\r\n% It is always a good idea to plot your data. In fact, I'll suggest\r\n% that you should plot everything. Plots are useful, since your eye\r\n% and brain are splendid at things like pattern recognition. Only\r\n% you know your data, as the scientist, engineer, or analyst. You\r\n% will always benefit if you can employ your knowledge of a system\r\n% as part of the modeling process.\r\n% \r\n% Next, always think in advance about your goals for any model.\r\n%\r\n% * Will you use this model purely for interpolation, i.e., for\r\n% predictive purposes only?\r\n%\r\n% * Do you need to derive some\r\n% understanding about the process from your model? Perhaps you need\r\n% to estimate an upper asymptote of your process. If so, then you\r\n% may want a model that has such an upper asymptote built into it.\r\n%\r\n% * Will you need to include the model coefficients in some paper that\r\n% you expect to write? Will you need to use this model in MATLAB, or\r\n% in some other tool?\r\n%\r\n% * Must the model be simple to evaluate?\r\n%\r\n% * Must the model be efficient to evaluate? You should know whether you will\r\n% need to evaluate this model millions\r\n% of times, or just once.\r\n%\r\n% * What do you know about your data? Is there\r\n% noise in the data? May you ignore that noise, or must you smooth\r\n% the noise away?\r\n%\r\n% * What do you know about the underlying functional relationship? Is\r\n% it monotone? Increasing? Decreasing? Positively\/negatively curved?\r\n% Must it pass through a specific point?\r\n%\r\n% * Must the interpolant have specific requirements in terms of\r\n% continuity or differentiability?\r\n% \r\n% For example, I once had a problem where I knew I had some\r\n% significant noise in our process, but I chose not to do any\r\n% smoothing anyway. Any such smoothing would also have smoothed\r\n% out some potentially important features of our process. Since I\r\n% could survive with the noise in my interpolant, I chose the lesser\r\n% evil.\r\n% \r\n% Always know your goals for any such task.\r\n\r\nplot(x,y,'bo')\r\nxlabel X\r\nylabel Y\r\ngrid on\r\ntitle 'Exponential data'\r\n\r\n%%\r\n% This is a nice, well-behaved function. It is of the form y = f(x). I'll\r\n% pretend for the moment that I have no idea what was in the underlying\r\n% functional relationship.\r\n%\r\n% One thing I learned in some long past calculus course was that a\r\n% Taylor series will provide an approximation for many functions.\r\n% Polynomials are useful things. They are simple to use, simple to\r\n% build, simple to work with. And a truncated Taylor series is\r\n% basically a nice, simple polynomial.\r\n% \r\n% So we will start with a linear polynomial approximation for this\r\n% curve, built using\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/polyfit.html |polyfit|>\r\n% . This is a utility provided in MATLAB to\r\n% estimate a polynomial model using linear regression techniques.\r\n% We could also use many other tools to build our polynomial model,\r\n% but polyfit is a useful one, and easy to use.\r\n% \r\n% A linear, or first\r\n% <http:\/\/en.wikipedia.org\/wiki\/Degree_of_a_polynomial degree>\r\n% polynomial (many use the words \"order\" and \"degree\" interchangeably),\r\n% might be written mathematically \r\n% as |y(x) = a1*x + a2|. In MATLAB we will merely store the coefficients,\r\n% as a vector [a1,a0]. Note that a polynomial in MATLAB has it's\r\n% coefficients stored with the highest order term first.\r\nP1 = polyfit(x,y,1)\r\n\r\n%%\r\n% We can evaluate the polynomial with \r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/polyval.html |polyval|>.\r\nyhat = polyval(P1,x);\r\n\r\nplot(x,y,'bo',x,yhat,'r-')\r\nxlabel X\r\nylabel Y\r\ngrid on\r\ntitle 'Linear polynomial fit'\r\n\r\n%% Look at the residuals\r\n% In fact, I'll claim the relationship we are modeling is not terribly well represented\r\n% by a linear model. Depending on your needs for this model, you might have\r\n% decided differently.\r\n%\r\n% When you build a regression model, look at the residuals.\r\n%\r\n%  ALWAYS PLOT EVERYTHING!\r\n%\r\n% Plot your residuals. In general, some good ways to plot the residuals are\r\n% versus\r\n%\r\n% * the independent variable. Look for patterns. Patterns here in the\r\n% residuals are often a hint that you should look more deeply.\r\n%\r\n% * the measurement order, just in case there was a problem with your\r\n% equipment. (I've seen many cases where measurement apparatus was\r\n% re-calibrated at the end of every day, but an experiment spanned more than\r\n% one day.)\r\n%\r\n% * the dependent variable. This might\r\n% help pick out cases of non-uniform variance.\r\n% \r\n% Look at your plots. Think about what you see there. Compare it to your expectations.\r\n% \r\n% Since this data was very simply generated, I'll dispense with some of\r\n% those plots for brevity.\r\n% Note that the residuals for this linear fit look vaguely like a\r\n% quadratic polynomial. This is often the case when there is lack\r\n% of fit in a polynomial. That lack of fit often looks like the first\r\n% term we truncated from the Taylor series.\r\nres = y - yhat;\r\nplot(x,res,'bo')\r\nxlabel X\r\nylabel Residuals\r\ngrid on\r\ntitle 'Residuals for the linear fit'\r\n\r\n%% \r\n% If the residuals looked vaguely parabolic in shape, then it might make sense to use\r\n% a second order (quadratic) polynomial for our fit.\r\nP2 = polyfit(x,y,2)\r\nyhat = polyval(P2,x);\r\nplot(x,y,'bo',x,yhat,'r-')\r\nxlabel X\r\nylabel Y\r\ngrid on\r\ntitle 'Quadratic polynomial fit'\r\n\r\n%%\r\n% Note that the residuals here look vaguely like a cubic polynomial,\r\n% although they are much smaller in magnitude than the previous fit. Again, at each step as\r\n% we increase the order of the model, the residuals will often tend to\r\n% look much like a polynomial of the next higher order.\r\nres = y - yhat;\r\nplot(x,res,'bo')\r\nxlabel X\r\nylabel Residuals\r\ngrid on\r\ntitle 'Residuals for the quadratic fit'\r\n\r\n%% Higher order polynomials - are more terms always better?\r\n% Polynomial modeling with polyfit is indeed simple and easy to do.\r\n% In fact, we might decide to just use higher and higher order polynomials,\r\n% always chasing the term we truncated in the previous model. After all,\r\n% if a quadratic model is better than a linear one, then why not go to a\r\n% cubic? Ten terms must be better than nine. Look at a tenth order model.\r\nP10 = polyfit(x,y,10);\r\ndisp(P10(1:6))\r\ndisp(P10(7:11))\r\n\r\n%%\r\n% We can write the\r\n% <http:\/\/en.wikipedia.org\/wiki\/Taylor_series |Maclaurin series|>\r\n% representation for the exponential function as\r\n% \r\n% $$exp(x) = 1 + x + x^2\/2 + x^3\/6 + x^4\/24 + x^5\/120 + x^6\/720 + x^7\/5040 + ...$$\r\n%\r\n% We can compare P10 to the coefficients of the known series. How well did\r\n% we do in the fit?\r\nseries = 1.\/factorial(10:-1:0);\r\ndisp(series(1:6))\r\ndisp(series(7:11))\r\n\r\n%%\r\n% Note that the higher order coefficients deviate somewhat from the known\r\n% series, although the lower order terms appear to be quite accurate.\r\nyhat = polyval(P10,x);\r\nplot(x,y,'bo',x,yhat,'r-')\r\nxlabel X\r\nylabel Y\r\ngrid on\r\ntitle 'Tenth order polynomial fit'\r\n\r\n%%\r\n% The residuals oscillate tightly around zero. In fact, they are so small\r\n% that this last polynomial begins to approach a true interpolating polynomial.\r\n% Perhaps if we just add a few more terms we may get there.\r\n% The numerical issues of floating point arithmetic will often preclude\r\n% true interpolation down to the least significant bit anyway.\r\nres = y - yhat;\r\nplot(x,res,'bo')\r\nxlabel X\r\nylabel Residuals\r\ngrid on\r\ntitle 'Residuals for the tenth order fit'\r\n\r\n%% What do you do with interpolation?\r\n% I'll start talking about true interpolation in my next blog. But remember that\r\n% interpolation is different from the approximations provided by polyfit or\r\n% any other regression modeling tool.\r\n% \r\n% Until that time, please give me your \r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=141#respond |comments|>\r\n% on this blog, or\r\n% ideas for future blog topics on interpolation or modeling in general.\r\n% Do you have some interesting applications of interpolation? Some\r\n% interesting problems?\r\n\r\n##### SOURCE END ##### 90e4d2c2328643a19d71a40dc659db85\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      I'd like to introduce a new guest blogger - John D'Errico - an applied mathematician, now retired from Eastman Kodak, where he used MATLAB for over 20 years. Since then, MATLAB is\r\n     ... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/11\/interpolation-in-matlab\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[23],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/141"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=141"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/141\/revisions"}],"predecessor-version":[{"id":2536,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/141\/revisions\/2536"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=141"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=141"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=141"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}