Turn off some warnings:

Off[LinearSolve::luc];

A Model Selection Task

You are given some data {((Subscript[x, 1]),(Subscript[y, 1])), ..., ((Subscript[x, n]),(Subscript[y, n]))} and you want to learn the mapping:

x ->y

or in particular the function g:

y = g(x)

for simplicity we will assume (Subscript[x, i]),(Subscript[y, i])∈ℝ.

Linear Models

Restrict the set of functions we can learn to:

(Subscript[f, w])(x)=(∑_(j=1)^d)(Subscript[w, j]) (Subscript[ϕ, j])(x)

Unobserved Process

Suppose nature creates the following function (which in practice we never get to see):

g(x) = 1 -2 x - 0.2 (x^2) + 0.2 (x^3)

Plotted as:

g = Function[x, 1 -2x - .2x^2 + 0.2 *x^3];
Plot[g[x],{x,-5,5}]

Synthetic Data

Instead we get to see noisy observations of that function:

(Subscript[y, i])= g((Subscript[x, i])) + (Subscript[ϵ, i])

The noise model is IMPORTANT:

(Subscript[ϵ, i])~N(0, (σ^2))

|(Subscript[f, w])(x) - y|

Plot[Exp[-Abs[x - 0]],{x,-5,5}]

for some unknown (σ^2) corresponds to standard Least squares regression:

(g((Subscript[x, i])) - (Subscript[y, i])) = (Subscript[ϵ, i])

n = 20;
data = Map[Function[x,{x,g[x] +Random[NormalDistribution[0,2]]}], RandomReal[{-5,5},n]];
ListPlot[data]

Basic Linear Regression

Lets try fitting a 10th degree polynomial:

d = 11;
features = Function[x,Table[(x^i),{i,0,d-1}]];
features[x]
Φ = Map[features,data[[All,1]]];
y = {data[[All,2]]}\[Transpose];

Solve for the weights

Algorithmic Note:

Solving is often more numerically stable than direct inversion:

(A^(-1))y = (Underscript[Solve[A,y], (Underscript[⏟, more stable])])

Φᵀ Φ w = Φᵀ y

w = (((Φᵀ Φ))^(-1))Φᵀ y

w = LinearSolve[Φ\[Transpose] .  Φ , Φ\[Transpose] .  y ];

fhat = Function[{x,w}, First[w\[Transpose] . features[x]]];
Show[Plot[fhat[x,w],{x,-5,5}],ListPlot[data]]

w

With regularization:

w = (((Φᵀ Φ + λ I))^(-1))Φᵀ y

Manipulate[
w = LinearSolve[(Φ\[Transpose] .  Φ + λ  IdentityMatrix[d]), Φ\[Transpose] .  y ];
{MatrixForm[w],
Show[ListPlot[data],Plot[{g[x],fhat[x,w]},{x,-5,5}]]},
{λ,0, 100},
SaveDefinitions->True
]

Different prior assumption:

Do we really believe that the intercept should have a zero mean prior:

Manipulate[
sigmaPrior = DiagonalMatrix[Join[{0},ConstantArray[1,d-1]]];
w = LinearSolve[(Φ\[Transpose] .  Φ + λ  sigmaPrior), Φ\[Transpose] .  y ];
{MatrixForm[w],
Show[ListPlot[data],Plot[{g[x],fhat[x,w]},{x,-5,5}]]},
{λ,0, 100},
SaveDefinitions->True
]

Infinite variance in the prior:

MatrixForm[sigmaPrior]

Thinking more about Features

We can consider other basis functions:

Manipulate[
Module[{Φ,w,f,sigmaPrior,ϕ},
ϕ = Join[{Function[x,1]},
Table[Function[x, Exp[((-(((x-b))^2))/σ)] ]/.b->i,{i,-5,5}]];
f = Function[x,(∑_(i=1)^(Length[ϕ])) ϕ[[i]][x]];
Plot[f[x],{x,-5,5}]
],
{σ,.01, 3}
]

Manipulate[
Module[{Φ,w,f,sigmaPrior,ϕ},
ϕ = Join[{Function[x,1]},
Table[Function[x, Exp[((-(((x-b))^2))/σ)] ]/.b->i,{i,-5,5}]];
sigmaPrior = DiagonalMatrix[Join[{0},ConstantArray[1,Length[ϕ]-1]]];
Φ= Map[Function[x,Map[Function[h,h[x]],ϕ]],data[[All,1]]];
w = LinearSolve[(Φ\[Transpose] .  Φ + λ  sigmaPrior), Φ\[Transpose] .  y ];
f =Function[x, (∑_(i=1)^(Length[ϕ]))w[[i]]ϕ[[i]][x]];
Show[ListPlot[data],Plot[{g[x],f[x]},{x,-5,5}]]
],
{λ,0, 10},
{σ,.01, 3}
]