Voorbeeld 1: lineaire regressie Eerst rekenen we het voorbeeld uit de lestekst na:

>>> import numpy as np
>>> import numpy.linalg as la
>>> A = np.array([[1,1],[1,2],[1,3]]); print(A) # matrix A voor regressielijn
[[1 1]
 [1 2]
 [1 3]]
>>> b = np.array([2,2,4]);  print(b)            # gegevensvector
[2 2 4]
>>> Ainv = la.pinv(A); print(np.round(Ainv, 4)) # pseudoinverse
[[ 1.3333  0.3333 -0.6667]
 [-0.5     0.      0.5   ]]
>>> print(np.round(Ainv @ b, 4))                # coëfficiënten regressielijn
[0.6667 1. ]     

Hier staat de uitkomst: \(y= 0.6667+ t\).
We kunnen het resultaat ook via de singulierewaardenontbinding berekenen.

>>> U, S, Vt = la.svd(A, full_matrices=False)   # SVD van A
>>> print(U)                  # matrix U in SVD
[[-0.32310986 0.85377594]
[-0.54750702  0.18321952]
[-0.77190419 -0.4873369 ]]
>>> print(np.diag(S))         # diagonaalmatrix met singuliere waarden
[[4.07914333 0.        ]
 [0.         0.60049122]]
>>> print(Vt)                 # getransponeerde van matrix V in SVD
[[-0.40266324 -0.91534819]
[ 0.91534819 -0.40266324]]
>>> print(U @ np.diag(S) @ Vt) # check van SVD berekening
[[1., 1.],
 [1., 2.],
 [1., 3.]])
>>> Splus = la.inv(np.diag(S)); print(Splus)
[[0.24514951, 0.        ]
 [0. ,        1.66530329]
>>> Ainv = Vt.transpose() @ Splus @ U.transpose()
>>> print(np.round(Ainv, 4))
[[ 1.3333 0.3333 -0.6667]
 [-0.5    0.     0.5    ]]
>>> print(np.round(Ainv @ b, 4))               # coëfficiënten regressielijn
[0.6667 1. ]     

Voorbeeld 2: kwadratische regressie Het tweede voorbeeld betreft een gegevensverzameling die we construeren door ruis aan de parabool \(1-2x+\frac{1}{2}x^2\) toe te voegen. We zoeken daarna de best bij de gegevens passende parabool. In onze dataset wordt deze bij benadering \(y=0.98-1.94x+0.49x^2\).

>>> import numpy as np
>>> import numpy.linalg as la
>>> import matplotlib.pyplot as plt
>>> x = np.arange(0, 4, 0.5)
>>> y = 1 - 2*x +1/2*x**2
>>> y = y + 0.1*max(y)*(2*np.random.randint(2, size=m) - np.ones(m))
>>> A = np.vstack((np.ones(m), x, x**2)).transpose()
>>> pinvA = la.pinv(A)
>>> c = pinvA @ y; print(c)
[ 0.975 -1.93571429 0.48571429]
>>> xx = np.linspace(-0.1,4.1,50)
>>> yy = c[0] + c[1]*xx + c[2]*xx**2
>>> plt.plot(x, y, 'bo')
>>> plt.plot(xx, yy, 'g-')
>>> plt.axis([-0.1, 4.1, -1.1, 1.1])
>>> plt.show()