Nichtlineare Regression, Klassifikation und Machine Learning - Übung 1

From StatWiki
Jump to: navigation, search

Crane, Protractor and River

Gegeben ist ein rechtwinkeliges Dreieck mit

  • Hypothenuse ist 21m lang
  • Die Messung der Länge der senkrechten Kante ist y_2 = 14.94 \mbox{m} \,
  • Die Messung des Winkels y_1 = 0.47 \, \mbox{rad} = \theta\,

Gesucht ist

  • Model
  • Compute linear approximation
  • Results of linear approximation for various \theta^{(0)} \, (Hinweis: y_1, y_2, \ldots \, sind gute Startpunkte)
  • Picture of expectation curve \{ \eta(\theta): \theta \in (0, \pi/2) \}\,
  • Formula and picture of K_{\mathrm{int}}, K_{\mathrm{par}} \,
  • Autofahrt überlegen
  • Bonus: \hat \theta \, finden (Hinweis: Minimum von S(\theta)\, unter Verwendung von numerischen Methoden; z.b. Bisektion)

Alternative Lösungen

Model

y = \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = \begin{pmatrix} \theta \\ 21 \, \sin(\theta) \end{pmatrix} = \begin{pmatrix} 0.47 \\ 14.94 \end{pmatrix}\,

y  = \eta(\theta) + \epsilon \,

\eta(\theta) = \begin{pmatrix} \theta \\ 21 \, \sin(\theta) \end{pmatrix}

Assertion: \operatorname{Var}(\epsilon) = \sigma^2 \,

Linear Approximation

\begin{matrix}
\underbrace{y - \eta(\theta^{(0)}) + \frac{d \, \eta(\theta^{(0)})}{d \, \theta} \eta(\theta^{(0)})} & = & \underbrace{\frac{d \, \eta(\theta^{(0)})}{d \, \theta}} & \theta + \epsilon \\
y^{\star} & & F & 
\end{matrix}

Linear approximation of \hat \theta\,: \hat{\theta}^{\star} = (F'F)^{-1} F'y^{\star} \,

F = \frac{d \, \eta(\theta^{(0)})}{d \, \theta} = \begin{pmatrix} 1 \\ 21 \ \cos(\theta^{(0)}) \end{pmatrix} \,

(F^TF)^{-1} = (1^2 + 21^2 \ \cos^2 (\theta^{(0)}))^{-1} (it's a number)

Results of linear approximation for various starting points

Measurement 1 in Degree: y_1 = 26.92902^{\circ}

Linear approximations

\begin{align}
\theta^{(0)} = y_1 & \Rightarrow \theta^{\star} = 18.73681^{\circ} \\
\theta^{(0)} = \frac{\arcsin(y_2)}{21} & \Rightarrow \theta^{\star} = 12.50937^{\circ} \\
\theta^{(0)} = 45^{\circ}& \Rightarrow \theta^{\star} = 12.50729^{\circ}
\end{align}

> y = c(.47, 14.94)
>
> options(warn = -1)
> library(circular)

Attaching package: ‘circular’

The following object is masked from ‘package:stats’:

    sd, var

>
> lin.approx = function(theta0) {
+         eta = function(theta0) c(theta0, 21 * sin(theta0))
+         d.eta = function(theta0) c(1, 21 * cos(theta0))
+
+         y.star = y - eta(theta0) + d.eta(theta0) * theta0
+
+         # eta star in degree
+         c(deg((d.eta(theta0) %*% d.eta(theta0))^-1 * d.eta(theta0) %*% y.star))
+ }
>
> deg(0.47)
[1] 26.92902
> lin.approx(0.47)
[1] 43.49673
> lin.approx(asin(14.94/21))
[1] 45.26707
> lin.approx(rad(45))
[1] 45.26702

Picture of expectation curve

> pdf(rpdf,width=12)
>
> drawPiAxis = function() {
+         at=seq(0,2,length.out=9)
+
+         # needed for fractions
+         library(MASS)
+
+         # attr needed to get e.g. 1/2 instead of 0.5
+         axis(1, at=at*pi, labels=sapply(attr(fractions(at),"fracs"), function(x) parse(text=paste(x, "* pi", sep=""))))
+ }
>
> eta=function(theta) c(theta, 21 * sin(theta))
>
> plot(t(sapply(seq(0,2*pi,length.out=1000), eta)),
+         main="Expectation Surface",
+         xlab=expression(y[1]), ylab=expression(y[2]),
+         type="l",
+         axes=F)
>
> drawPiAxis()
> axis(2, c(-21,-15,-10,-5,0,5,10,15,21))
>
> text(pi/8, 15, expression( eta(theta) : theta %in% (list(0, pi/2)) ))

Formula and picture of curvature

Intrinsic curvature K_{\mbox{int}}(\theta) := \frac{ || (I-P(\theta)) \frac{ \partial^2 \eta(\theta) } {\partial^2 \theta^T} || } { || \frac{\partial \eta(\theta) }{\partial \theta^T} ||^2 }

Parameter effect curvature K_{\mbox{par}}(\theta) := \frac{ || P(\theta) \frac{ \partial^2 \eta(\theta) } {\partial^2 \theta^T} || } { || \frac{\partial \eta(\theta) }{\partial \theta^T} ||^2 }

where P(\theta) = F(F^T F)^{-1} F^T\, the orthogonal project is.

Interpretation

K_{\mbox{int}}(\theta)\, ist die Inverse des Radius jenes Kreises, der sich an der Stelle \theta\, optimal an die Funktion \eta(\theta)\, anpasst. 
Je grösser der Radius des besser kann linear approxmiert werden.

K_{\mbox{par}}(\theta)\, gibt das Ausmass der Änderung von \eta(\theta)\, an, wenn \theta\, geändert wird.

Hier sind 2 Fälle zu unterscheiden

1. K_{\mbox{par}}(\theta) = 0\, heisst, dass die Abstände zwischen aufeinander folgenden Punkten |\eta(\theta_i) - \eta(\theta_{i+1})|\, konstant ist.
2. K_{\mbox{par}}(\theta) > 0\, heisst, dass die Abstände zwischen aufeinander folgenden Punkten |\eta(\theta_i) - \eta(\theta_{i+1})|\, unterschiedlich sind. Das ist der schwierigere Fall.

\eta(\theta) = \begin{pmatrix} \theta \\ 21 \, \sin(\theta) \end{pmatrix}

F =  \frac{d \, \eta(\theta)}{d \, \theta} = \begin{pmatrix} 1 \\ 21 \ \cos(\theta) \end{pmatrix} \,

(F^TF)^{-1} = \frac{1}{1 + (21 \ \cos (\theta))^2}

Projector

 \begin{align}
P(\theta)
& = F (F^TF)^{-1} F^T \\
& =  \begin{pmatrix} 1 \\ 21 \ \cos(\theta) \end{pmatrix}  \frac{1}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} 1 \\ 21 \ \cos(\theta) \end{pmatrix}^T \\
& =  \frac{1}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} 1 \\ 21 \ \cos(\theta) \end{pmatrix}  \begin{pmatrix} 1 \\ 21 \ \cos(\theta) \end{pmatrix}^T \\
& =  \frac{1}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} 1 & 21 \ \cos(\theta) \\ 21 \ \cos(\theta) & (21 \ \cos(\theta))^2  \end{pmatrix} \\
\end{align}


 \begin{align}
I - P(\theta)
& = \frac{1}{1 + (21 \ \cos(\theta))^2} \begin{pmatrix} 1 + (21 \ \cos(\theta))^2 & 0 \\ 0 & 1 + (21 \ \cos(\theta))^2 \end{pmatrix} - P \\
& = \frac{1}{1 + (21 \ \cos(\theta))^2} \begin{pmatrix} 1 + (21 \ \cos(\theta))^2 & 0 \\ 0 & 1 + (21 \ \cos(\theta))^2 \end{pmatrix} - \frac{1}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} 1 & 21 \ \cos(\theta) \\ 21 \ \cos(\theta) & (21 \ \cos(\theta))^2  \end{pmatrix} \\
& = \frac{1}{1 + (21 \ \cos(\theta))^2} \left( \begin{pmatrix} 1 + (21 \ \cos(\theta))^2 & 0 \\ 0 & 1 + (21 \ \cos(\theta))^2 \end{pmatrix} - \begin{pmatrix} 1 & 21 \ \cos(\theta) \\ 21 \ \cos(\theta) & (21 \ \cos(\theta))^2  \end{pmatrix} \right) \\
& = \frac{1}{1 + (21 \ \cos(\theta))^2} \begin{pmatrix} (21 \ \cos(\theta))^2 & -21 \ \cos(\theta) \\ -21 \ \cos(\theta) & 1 \end{pmatrix} \\
\end{align}
 \frac{d^2 \, \eta(\theta)}{d^2 \, \theta} = \begin{pmatrix} 0 \\ - 21 \ \sin(\theta) \end{pmatrix} \,
 \begin{align}
\left|\left| (I - P(\theta)) \frac{d^2 \, \eta(\theta)}{d^2 \, \theta} \right|\right|
& = \left|\left| \frac{1}{1 + (21 \ \cos(\theta))^2} \begin{pmatrix} (21 \ \cos(\theta))^2 & -21 \ \cos(\theta) \\ -21 \ \cos(\theta) & 1 \end{pmatrix} \begin{pmatrix} 0 \\ - 21 \ \sin(\theta) \end{pmatrix} \right|\right| \\ 
& = \left|\left| \frac{1}{1 + (21 \ \cos(\theta))^2} \begin{pmatrix} 21^2 \ \cos(\theta) \sin(\theta) \\ -21 \ \sin(\theta) \end{pmatrix} \right|\right| \\
& = \left|\left| \frac{-21 \ \sin(\theta) }{1 + (21 \ \cos(\theta))^2} \begin{pmatrix} 21 \ \cos(\theta) \\ 1 \end{pmatrix} \right|\right| \\
& = \left| \frac{-21 \ \sin(\theta) }{1 + (21 \ \cos(\theta))^2} \right| \sqrt{ 1^2 + (21 \ \cos(\theta))^2 } \\
& = \frac{21 \ |\sin(\theta)| }{\sqrt{1 + (21 \ \cos(\theta))^2}} \\
\end{align}
\begin{align}
K_{\mbox{int}}(\theta) 
& = \frac{ || (I-P(\theta)) \frac{ \partial^2 \eta(\theta) } {\partial^2 \theta^T} || } { || \frac{\partial \eta(\theta) }{\partial \theta^T} ||^2 } \\
& = \frac{ \frac{21 \ |\sin(\theta)| }{\sqrt{1 + (21 \ \cos(\theta))^2}} } { 1 + (21 \ \cos(\theta))^2 } \\
& = \frac{21 \ |\sin(\theta)| }{ (1 + (21 \ \cos(\theta))^2)^{\frac 1 2} (1 + (21 \ \cos(\theta))^2) } \\
& = \frac{21 \ |\sin(\theta)| }{ (1 + (21 \ \cos(\theta))^2)^{1 + \frac 1 2} } \\
& = \frac{21 \ |\sin(\theta)| }{ (1 + (21 \ \cos(\theta))^2)^{\frac 3 2} }
\end{align}
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/671756eba1f43acb36337c82d26069ff6ec6918b_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf,width=12)
> x=seq(0,2*pi,length.out=1000)
> kint = function(x) (21^2 * abs(sin(x)))/(1 + (21*cos(x))^2)^(3/2)
> plot(x, kint(x),type="l",axes=F,xlab=expression(theta),ylab=expression(K[int](theta)),main="Intrinsic Curvature")
>
> drawPiAxis()
Error: could not find function "drawPiAxis"
Execution halted
in

pdf(rpdf,width=12) x=seq(0,2*pi,length.out=1000) kint = function(x) (21^2 * abs(sin(x)))/(1 + (21*cos(x))^2)^(3/2) plot(x, kint(x),type="l",axes=F,xlab=expression(theta),ylab=expression(K[int](theta)),main="Intrinsic Curvature")

drawPiAxis() axis(2)

\begin{align}
\left|\left| P(\theta) \frac{d^2 \, \eta(\theta)}{d^2 \, \theta} \right|\right|
& = \left|\left| \frac{1}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} 1 & 21 \ \cos(\theta) \\ 21 \ \cos(\theta) & (21 \ \cos(\theta))^2  \end{pmatrix} \begin{pmatrix} 0 \\ - 21 \ \sin(\theta) \end{pmatrix} \right|\right| \\
& = \left|\left| \frac{1}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} - 21^2 \ \sin(\theta) \cos(\theta) \\ -21^3 \ \sin(\theta) \cos^2(\theta) \end{pmatrix} \right|\right| \\
& = \left|\left| \frac{- 21^2 \ \sin(\theta) \ \cos(\theta)}{1 + (21 \ \cos (\theta))^2} \begin{pmatrix} 1 \\ -21 \ \cos(\theta) \end{pmatrix} \right|\right| \\
& = \frac{21^2 |\sin(\theta)| |\cos(\theta)|}{1 + (21 \ \cos (\theta))^2} \sqrt{1^2 + (21 \ \cos(\theta))^2} \\
& = \frac{21^2 |\sin(\theta)| |\cos(\theta)|}{\sqrt{1 + (21 \ \cos(\theta))^2}}
\end{align}
\begin{align}
K_{\mbox{par}}(\theta) 
& = \frac{ || (I-P(\theta)) \frac{ \partial^2 \eta(\theta) } {\partial^2 \theta^T} || } { || \frac{\partial \eta(\theta) }{\partial \theta^T} ||^2 } \\
& = \frac{ \frac{21^2 |\sin(\theta)| |\cos(\theta)|}{\sqrt{1 + (21 \ \cos(\theta))^2}} }{ 1 + (21 \ \cos(\theta))^2 } \\
& = \frac{21^2 |\sin(\theta)| |\cos(\theta)|}{(1 + (21 \ \cos(\theta))^2)^{\frac 3 2}} \\
\end{align}
REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/cf4a243b55a0d940b4d8019f2c9fda3180f393fc_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf,width=12)
> x=seq(0,2*pi,length.out=1000)
> kpar = function(x) (21^2 * abs(sin(x)) * abs(cos(x)))/(1 + (21*cos(x))^2)^(3/2)
> plot(x, kpar(x),type="l",axes=F,xlab=expression(theta),ylab=expression(K[int](theta)),main="Parameter Effect Curvature")
>
> drawPiAxis()
Error: could not find function "drawPiAxis"
Execution halted
in

pdf(rpdf,width=12) x=seq(0,2*pi,length.out=1000) kpar = function(x) (21^2 * abs(sin(x)) * abs(cos(x)))/(1 + (21*cos(x))^2)^(3/2) plot(x, kpar(x),type="l",axes=F,xlab=expression(theta),ylab=expression(K[int](theta)),main="Parameter Effect Curvature")

drawPiAxis() axis(2)

Autofahrt überlegen

The parameter effect curvature represents a measure of change of the intrinsic curvature.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/40d6ea148a4a21d98c0ba7f85e9c722d499781e6_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf,width=12,height=12)
>
> par(mfrow=c(3,1))
> plot(t(sapply(seq(0,2*pi,length.out=1000), eta)),
+ 	main="Expectation Surface",
+ 	xlab=expression(y[1]), ylab=expression(y[2]),
+         type="l",
+         axes=F)
Error in match.fun(FUN) : object 'eta' not found
Calls: plot -> t -> sapply -> match.fun
Execution halted
in

pdf(rpdf,width=12,height=12)

par(mfrow=c(3,1)) plot(t(sapply(seq(0,2*pi,length.out=1000), eta)), main="Expectation Surface", xlab=expression(y[1]), ylab=expression(y[2]),

       type="l", 
       axes=F)

drawPiAxis() axis(2, c(-21,-15,-10,-5,0,5,10,15,21))

x=seq(0,2*pi,length.out=1000)

plot(x, kint(x),type="l",axes=F,xlab=expression(theta),ylab=expression(K[int](theta)),main="Intrinsic Curvature") drawPiAxis() axis(2)

plot(x, kpar(x),type="l",axes=F,xlab=expression(theta),ylab=expression(K[int](theta)),main="Parameter Effect Curvature") drawPiAxis() axis(2)


Minimum von S(.) finden (Bonus)

S(\theta) := || y - \eta(\theta) ||^2 \,

\hat \theta = \arg\min S(\theta) \,

Minimum durch  \frac{d S(\theta)}{\theta} = 0 \,

S(\theta) = (y_1 - x)^2 + (y_2 - 21 sin(x))^2 \,

<Axiom echo="true"> D((y[1] - x)^2 + (y[2] - 21 * sin(x))^2, x) </Axiom>

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/ab6beb046d2bc1d02a07166e98c5a8d3b945c00f_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf)
> s = function(x) (.47 - x)^2 + (14.94 - 21*sin(x))^2
> x=seq(0,pi/2,length.out=1000)
> plot(x, s(x),type="l",axes=F,xlab=expression(theta),ylab=expression(S(theta)),main="S")
> drawPiAxis()
Error: could not find function "drawPiAxis"
Execution halted
in

pdf(rpdf) s = function(x) (.47 - x)^2 + (14.94 - 21*sin(x))^2 x=seq(0,pi/2,length.out=1000) plot(x, s(x),type="l",axes=F,xlab=expression(theta),ylab=expression(S(theta)),main="S") drawPiAxis() axis(2)

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/bfab40598726c9e5f3873c10ff353d5dc4f00dda_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf)
> s.derived = function(x) 21*21 *cos(x) *sin(x) - 21 * 14.94 * cos(x) + x - .47
>
> x=seq(0,pi/2,length.out=1000)
> plot(x, s.derived(x),type="l",axes=F,xlab=expression(theta),ylab=expression(partialdiff * S(theta)/partialdiff * theta),main="Derivate of S")
> drawPiAxis()
Error: could not find function "drawPiAxis"
Execution halted
in

pdf(rpdf) s.derived = function(x) 21*21 *cos(x) *sin(x) - 21 * 14.94 * cos(x) + x - .47

x=seq(0,pi/2,length.out=1000) plot(x, s.derived(x),type="l",axes=F,xlab=expression(theta),ylab=expression(partialdiff * S(theta)/partialdiff * theta),main="Derivate of S") drawPiAxis() axis(2)

abline(h=0,col=2)

Find  \frac{d S(\theta)}{\theta} = 0 \, and get \hat \theta = 45.26273^\circ

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/6936d2265ff546baf53f4547dba6d6182245d292_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'text'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> options(warn=-1)
> library(circular)

Attaching package: ‘circular’

The following object is masked from ‘package:stats’:

    sd, var

>
> deg(uniroot(s.derived, c(0,pi/2))$root)
Error in uniroot(s.derived, c(0, pi/2)) : object 's.derived' not found
Calls: deg -> uniroot
Execution halted
in

options(warn=-1) library(circular)

deg(uniroot(s.derived, c(0,pi/2))$root)

Alternative \arg\min S(\theta) \, and get \hat \theta = 45.26724^\circ

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/5ea446e3655a62f59ec19f4e645776f0e4cda488_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'text'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> options(warn=-1)
> library(circular)

Attaching package: ‘circular’

The following object is masked from ‘package:stats’:

    sd, var

>
> deg(optimise(s, c(0,pi/2), maximum=F)$minimum)
Error in (function (arg)  : object 's' not found
Calls: deg -> optimise -> <Anonymous>
Execution halted
in

options(warn=-1) library(circular)

deg(optimise(s, c(0,pi/2), maximum=F)$minimum)