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

> library(pracma)
A <- matrix(c(1, 1, 1, 1, 2, 3), nrow=3); A   # matrix A voor regressielijn
     [,1] [,2]
[1,]    1    1
[2,]    1    2
[3,]    1    3
> b <- c(2, 2, 4)                  # gegevensvector
b =
     2
     2
     4
> Ainv <- round(pinv(A), 4); Ainv  # pseudoinverse
        [,1]   [,2]    [,3]
[1,]  1.3333 0.3333 -0.6667
[2,] -0.5000 0.0000  0.5000 
> Ainv %*% b                       # coëfficiënten regressielijn
       [,1]
[1,] 0.6667
[2,] 1.0000 

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

> svdA <- svd(A)  # singulierewaardenontbinding
> U <- svdA$u; U
           [,1]       [,2]
[1,] -0.3231099  0.8537759
[2,] -0.5475070  0.1832195
[3,] -0.7719042 -0.4873369
> S <- diag(svdA$d); S
         [,1]      [,2]
[1,] 4.079143 0.0000000
[2,] 0.000000 0.6004912
> V <- svdA$v; V
           [,1]       [,2]
[1,] -0.4026632  0.9153482
[2,] -0.9153482 -0.4026632
> Splus <- inv(S); Splus
          [,1]     [,2]
[1,] 0.2451495 0.000000
[2,] 0.0000000 1.665303
> Ainv <- V %*% Splus %*% t(U); Ainv  # pseudoinverse
          [,1]          [,2]       [,3]
[1,]  1.333333  3.333333e-01 -0.6666667
[2,] -0.500000 -4.163336e-17  0.5000000
> Ainv %*% b
          [,1]
[1,] 0.6666667
[2,] 1.0000000

Voorbeeld 2: kwadratische regressie Het tweede voorbeeld betreft een gegevensverzameling die we construëren 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.95-2.06x+0.52x^2\).

> library(pracma)
> x <- seq(from=0, to=4, by=0.5)
> m <- length(x)
> y <- rep(1,m) - 2*x + x^2/2   # parabool
> y <- y + 0.1*max(y)*sample(c(-1,1), m, replace=TRUE)  # 10 % ruis toegevoegd
> A = cbind(rep(1,m), x, x^2)
> pinvA <- pinv(A)
> c <- pinvA %*% y; c # pseudoinverse van A geeft regressiecoëfficiënten
           [,1]
[1,]  0.9545455
[2,] -2.0601732
[3,]  0.5233766
> xx <- seq(from=-0.1, to=4.1, length=50)
> yy <- rep(c[1],50) + c[2]*xx + c[3]*xx^2
> plot(x, y, type="p", col="blue", pch=1)
> lines(xx, yy, col="green", lwd=2) # visualisatie van data + benadering