# 1 Curve Fitting (NumPy)

## 1.1 Data generation
Let the following polynomial be the "true" model. 

\begin{equation}
t = w_2x^2 + w_1x +w_0
\end{equation}


First we sample $n$ equidistant values for $x$, in a range of $[0,1]$; Then compute $n$ corresponding target values for $t$ using the above polynomial equation; To "simulate" the real-world scenario where the observed data usually come with noise, add random Gaussian noise to our target values.


```python
import random


# get ground-truth data from the "true" model 
n = 20  # number of data samples
x = [(idx-round(n/2))/(n/2) for idx in range(n)]
print(x)

w = [1, 2, 3]
t = [xn**2*w[2] + xn*w[1] + w[0] for xn in x]
print(t)

# "simulated" observed target values
std_noise = 0.2
t_observed = [t[idx]+random.gauss(0,std_noise) for idx in range(n)]
```

    [-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
    [2.0, 1.6300000000000001, 1.3200000000000003, 1.0699999999999998, 0.8800000000000001, 0.75, 0.68, 0.67, 0.72, 0.83, 1.0, 1.23, 1.52, 1.87, 2.2800000000000002, 2.75, 3.2800000000000002, 3.8699999999999997, 4.5200000000000005, 5.23]


Python plotting library [matplotlib](https://matplotlib.org/) is used to plot the observed target values *t_observed* with respect at each sampling location $x$ in blue points [$x$,*t_observed*]. We can also plot noise-free "true" target values $t$ versus the same $x$. This represents the underlying model as a curve, hence *curve fitting*.


```python
# this line is useful for jupyter notebook only
%matplotlib inline

import matplotlib.pyplot as plt


# plot the curve and the noise-corrupted data
plt.plot(x,t,'r')
plt.plot(x,t_observed,'bo')
```

## 1.2 Model Fitting
For numerical computations such as curve fitting, we can use the powerful [numpy](http://www.numpy.org/). In this tutorial, we use [*numpy.linalg.lstsq*](https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsq) to solve [a system of linear equations](https://en.wikipedia.org/wiki/System_of_linear_equations) $\textbf{a}\textbf{x}=\textbf{b}$ for a [least sqaures solution](https://en.wikipedia.org/wiki/Least_squares).

Re-arrange the observed data to a linear system in matrix form "$\textbf{a}\textbf{x}=\textbf{b}$" (N.B. the $\textbf{x}$ is the unknown here, not the learning input $x$):

\begin{equation}
\textbf{t} = \textbf{X}\textbf{w}
\end{equation}

that is:

\begin{equation}
\begin{vmatrix}
t_{(1)} \\ t_{(2)} \\ \cdots \\ t_{(n)}
\end{vmatrix} = \begin{vmatrix}
x^2_{(1)} & x_{(1)} & 1 \\
x^2_{(2)} & x_{(2)} & 1 \\
& \cdots \\
x^2_{(n)} & x_{(n)} & 1 
\end{vmatrix} \cdot
\begin{vmatrix}
w_2 \\ w_1 \\ w_0
\end{vmatrix}
\end{equation}




```python
import numpy as np


# then use nunmpy for a least-square solution to the linear system "Xw=t"
t_observed = np.reshape(t_observed, [-1, 1])
x_1 = np.reshape(x, [-1, 1])
x_2 = np.square(x_1)
x_0 = np.ones_like(x_1)
X = np.concatenate([x_2,x_1,x_0],1)
# print to check the inputs
print(X)
print(X.shape)
print(t_observed.shape)
w_estimate = np.linalg.lstsq(X, t_observed, rcond=None)
print(w_estimate[0])  # print the output

# plot to see the estimated curve, i.e.
# t_estimate = [xn**2*w_estimate[2]+xn*w_estimate[1]+w_estimate[0] for xn in x]
# but matrix multiplication is more compact:
t_estimate = np.matmul(X,w_estimate[0])
```

    [[ 1.   -1.    1.  ]
     [ 0.81 -0.9   1.  ]
     [ 0.64 -0.8   1.  ]
     [ 0.49 -0.7   1.  ]
     [ 0.36 -0.6   1.  ]
     [ 0.25 -0.5   1.  ]
     [ 0.16 -0.4   1.  ]
     [ 0.09 -0.3   1.  ]
     [ 0.04 -0.2   1.  ]
     [ 0.01 -0.1   1.  ]
     [ 0.    0.    1.  ]
     [ 0.01  0.1   1.  ]
     [ 0.04  0.2   1.  ]
     [ 0.09  0.3   1.  ]
     [ 0.16  0.4   1.  ]
     [ 0.25  0.5   1.  ]
     [ 0.36  0.6   1.  ]
     [ 0.49  0.7   1.  ]
     [ 0.64  0.8   1.  ]
     [ 0.81  0.9   1.  ]]
    (20, 3)
    (20, 1)
    [[3.04570644]
     [1.85739825]
     [0.97774253]]


*w_estimate* are the best estimates of the weights $w$ - in a least squares sense. Finally, we can add a green curve (formed by points [$x$, *t_estimate*]) to the previous plot to visualise the estimated polynomial.


```python
plt.plot(x,t,'r')
plt.plot(x,t_observed,'bo')
plt.plot(x,t_estimate,'g')
```

## 1.3 Model Fitting Error
From the above plot, we can see that the green curve (i.e. the fitted polynomial model) does not go through all the observed points (in blue) which were used in fitting the model. This is known as redidual error, also known as training error. There is also discrepancy between the "true" curve (in red) and the estimated one (in green), but this difference is often impossible to obtain without access to the true underlying model.


```python
# residuals:
Residuals = t_estimate-t_observed
SR = np.sum(np.square(Residuals))  # sums of residuals: b - a*x
# root-mean-square error
RMSE = np.sqrt(np.mean(np.square(Residuals)))
print(SR)
print(RMSE)

# plot the error distribution
plt.hist(Residuals)
```

## Questions
### Sample size
- Effect on changing the sample size on model fitting and its errors.
- How many samples are needed? 

### Model fitting
- What is the objective (loss) function in this curve fitting problem?
- How does change of the noise level (i.e. *std_noise*) affect the model fitting?
- Are residualds a good estimate of the difference to the "true" target values?
- How to measure how good the estimates are?

852f6a215c6e09a84da99bae672499fc2cdb26c3
ipynb
Jupyter Notebook
tutorials/regression/01-CurveFitting.ipynb
YipengHu/MPHY0043
80aee58801def8fd1176fbc1a1f1892a9a98e172
Apache-2.0
tutorials/regression/01-CurveFitting.ipynb
YipengHu/MPHY0043
80aee58801def8fd1176fbc1a1f1892a9a98e172
Apache-2.0
tutorials/regression/01-CurveFitting.ipynb
YipengHu/MPHY0043
80aee58801def8fd1176fbc1a1f1892a9a98e172
Apache-2.0
2021-12-15T12:22:24.000Z
2021-12-15T12:22:24.000Z
Qwen/Qwen-72B
1. YES
2. YES