{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"convergence achieved in 5 iterations\n",
" initial final +/-\n",
" 0.500000 1.169004 0.047758\n",
" 1.000000 1.898204 0.065892\n"
]
}
],
"source": [
"% a minimal test of least squares fitting in octave\n",
"\n",
"SNR=100/3; % 3% noise\n",
"a=1.2; % \"true\" values\n",
"b=1.8;\n",
"\n",
"x=linspace(0, 5); \n",
"y = 1./(1 + a*x.^b) + randn(1,100)*0.05;\n",
"\n",
"pkg load optim\n",
"% fit using leasqr (x, y, pin, F), where\n",
"% x - Vector (or matrix) of independent variables.\n",
"% y - Vector or matrix of observed values.\n",
"% pin - Vector of initial parameters to be adjusted by leasqr.\n",
"% F - Name of function or function handle. \n",
"% The function must be of the form y = f(x, p), with y, x, p of the form y, x, pin.\n",
"% Outputs: [f, p, cvg, iter, corp, covp, covr, stdresid, Z, r2] = leasqr (…)\n",
"% f - Column vector of values computed: f = F(x,p).\n",
"% p - Column vector trial or final parameters, i.e, the solution.\n",
"% cvg - Scalar: = 1 if convergence, = 0 otherwise.\n",
"% iter - Scalar number of iterations used.\n",
"% corp - Correlation matrix for parameters.\n",
"% covp - Covariance matrix of the parameters.\n",
"% ...\n",
"% A simple-minded interpretation of the results of the fit:\n",
"% From http://originlab.com/www/helponline/Origin/en/UserGuide/The_Fit_Results.html\n",
"% \"The square root of a main diagonal value of the covariance matrix is the \n",
"% Standard Error of the corresponding parameter\"\n",
"\n",
"function y = f(x, p)\n",
" y = 1./(1+p(1)*x.^p(2));\n",
"end\n",
"\n",
"p_in = [0.5;1];\n",
"[yfit p_out cvg iter corp covp] = leasqr(x, y, p_in, \"f\");\n",
"\n",
"if (cvg == 1) \n",
" disp(sprintf(\"convergence achieved in %d iterations\\n initial final +/-\",iter))\n",
" er=sqrt(diag(covp));\n",
" disp([p_in p_out er])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%plot --format svg -w 600 -h 400"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
" figure(1); clf(1);\n",
" hold on;\n",
" plot(x,y,...\n",
" sprintf(\";input data, true a=%.2f, b=%.2f;b-\",a,b));\n",
" plot(x,f(x,p_in),...\n",
" sprintf(\";initial guess: a=%.2f, b=%.2f;g-\",p_in(1),p_in(2)));\n",
" plot(x,f(x,p_out),...\n",
" sprintf(\";best fit: a=%.2f +/- %.2f, b=%.2f +/- %.2f;r-\",p_out(1),er(1),p_out(2),er(2)));\n",
" dy=y-f(x,p_out);\n",
" plot(x,dy,...\n",
" sprintf(\";misfit, total misfit^2=%.4f;c-\",sum(dy.^2))); % could also use sumsq(dy)\n",
"\n",
" xlabel(\"x\");\n",
" ylabel(\"y = [ 1 + ax^b ]^{-1}\");\n",
" hold off;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Octave",
"language": "octave",
"name": "octave"
},
"language_info": {
"file_extension": ".m",
"help_links": [
{
"text": "GNU Octave",
"url": "https://www.gnu.org/software/octave/support.html"
},
{
"text": "Octave Kernel",
"url": "https://github.com/Calysto/octave_kernel"
},
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "octave",
"version": "4.2.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}