Untitled

1118 days ago by pub

Example 2:

A function of the form $f(x,{a_1},{a_2}) = \Large{\frac{a_1}{1+{a_2}e^{-x} } }$  yields the given measured data:

x=(-1, -0.6, -0.2, 0.2, 0.6, 1, 1.4, 1.8, 2.2, 2.6, 3)
y=(0.4, 0.53, 0.72, 0.89, 1.11, 1.23, 1.53, 1.6, 1.7, 1.85, 1.87).

Use the method of least squares to find the constants $a_1$ and $a_2$ so that the function best matches the data.

Find the standard deviation of the resulting funtion and the given data.

Solution:

Using the statistical definition of n, there are 11 pieces of data so n=11.

The number of constants is 2, that is, $a_1$ and $a_2$ so m=2

We then calculate the standard deviation using the statistical definition of n=number of points, m=number of constants: $\sigma=\large{\sqrt{\frac{S}{n-m}}}$ where $S=\sum (y2-y1)^2$.

In order to create our matrix equation, we must somehow make f(x) into a function with "separate" constants. Keep reading and you will see what we mean by this.

We look carefully at the function and see that we can make "separate" constants by inverting the function.

$\large{ \frac{1}{y} = \frac{1}{a_1} \cdot (1+a_2 e^{-x}) = \frac{1}{a_1} +\frac{a_2}{a_1} e^{-x} }$

We see that we can put this in the form $a+b e^{-x}$  where $a=\large{ \frac{1}{a_1} }$ and $b=\large{\frac{a_2}{a_1} }$ and .

This means that when we find a and b, we calculate: ${a_1}=1/a$ and ${a_2}=b \cdot a_1=b/a$.

We need to solve the (non-square) matrix equation AX=B using the least squares method where

 A is the nx2 dimensional array $A=\left[ 1, e^{-x} \right]$ X is the 2x1 dimensional array $\left[ {\begin{array}{*{20}{c} }a \\ b \end{array} } \right]$ and B is the nx1 dimension array $B=\large{\frac{1}{y}}$ AX=B $\left[ {\begin{array}{*{20}{c}} {1}&e^{1}\\ {1}&e^{0.6}\\ {1}&e^{0.2}\\ {1}&e^{-0.2} \\{1}&e^{-0.6} \\ {1}&e^{-1} \\ {1}&e^{-1.4} \\ {1}&e^{-1.8} \\ {1}&e^{-2.2} \\ {1}&e^{-2.6} \\ {1}&e^{-3} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c} }a \\ b \end{array} } \right] = \left[ {\begin{array}{*{20}{c}}{0.4} \\{0.53} \\{0.72} \\{0.89} \\{1.11} \\{1.23} \\{1.53} \\{1.6} \\{1.7} \\{1.85} \\ {1.87} \end{array}} \right]$
# In the applet below, the sliders are set at values approximating the "good values", that is, where the standard deviation is minimum. Click and drag them to other values and then see if you can get them "right" again. (If you look at the applet itself, n=number of subintervals and so in the definition of sigma (E13), we use n+1-m)

import numpy as np import numpy.linalg

The syntax for x1 and y1 are easy since they just need to be 1xn arrays. (We save x and y for variable names for plots.)

x1=np.array([-1, -0.6, -0.2, 0.2, 0.6, 1, 1.4, 1.8, 2.2, 2.6, 3]); y1=np.array([0.4, 0.53, 0.72, 0.89, 1.11, 1.23, 1.53, 1.6, 1.7, 1.85, 1.87]); print x1 print y1
 [-1. -0.6 -0.2 0.2 0.6 1. 1.4 1.8 2.2 2.6 3. ] [ 0.4 0.53 0.72 0.89 1.11 1.23 1.53 1.6 1.7 1.85 1.87] [-1. -0.6 -0.2 0.2 0.6 1. 1.4 1.8 2.2 2.6 3. ] [ 0.4 0.53 0.72 0.89 1.11 1.23 1.53 1.6 1.7 1.85 1.87]

The syntax for the array A is VERY important. It must be an nx2 array.

n=len(x1) A=np.array([[1, exp(-x1[j])] for j in range(n)]) print A
 [[ 1. 2.71828183] [ 1. 1.8221188 ] [ 1. 1.22140276] [ 1. 0.81873075] [ 1. 0.54881164] [ 1. 0.36787944] [ 1. 0.24659696] [ 1. 0.16529889] [ 1. 0.11080316] [ 1. 0.07427358] [ 1. 0.04978707]] [[ 1. 2.71828183] [ 1. 1.8221188 ] [ 1. 1.22140276] [ 1. 0.81873075] [ 1. 0.54881164] [ 1. 0.36787944] [ 1. 0.24659696] [ 1. 0.16529889] [ 1. 0.11080316] [ 1. 0.07427358] [ 1. 0.04978707]]

B is an nx1 array, but as in our first example, we will write it as a 1xn array. This yields X as a 1xn array.

B=1/y1 print B
 [ 2.5 1.88679245 1.38888889 1.12359551 0.9009009 0.81300813 0.65359477 0.625 0.58823529 0.54054054 0.53475936] [ 2.5 1.88679245 1.38888889 1.12359551 0.9009009 0.81300813 0.65359477 0.625 0.58823529 0.54054054 0.53475936]

Now we want to solve the (non-square) matrix equation AX=B using least squares. In MatLab we could just write X=A\B.

Here we use the numpy linalg command lstsq. Reference (Remember: X is just the constants $a$ and $b$.)

Notice the syntax that yields a quick way to get to the constants a and b here.

a,b=np.linalg.lstsq(A,B) print a,b
 0.500749759314 0.742519612115 0.500749759314 0.742519612115
a1=1/a a2=b/a print a1, a2
 1.99700545312 1.48281571444 1.99700545312 1.48281571444
var ( 'x' ) P=list_plot(zip(x1,y1)) C=plot(a1/(1+a2*exp(-x)),(x,-4,4), color='red', figsize=4) show(P+C)  y2=a1/(1+a2*exp(-x1)) print y2
 [ 0.39696286 0.53945908 0.71039616 0.90197889 1.10101461 1.29214415 1.46230291 1.60388158 1.71519738 1.7988868 1.85971212] [ 0.39696286 0.53945908 0.71039616 0.90197889 1.10101461 1.29214415 1.46230291 1.60388158 1.71519738 1.7988868 1.85971212]
S=sum((y2-y1)^2) print S
 0.0118243818139 0.0118243818139
m=2 sigma=sqrt(S/(n-m)) print sigma
 0.0362466577983 0.0362466577983