Problem 1

(1)

# The function f(x1, x2) = cos(x1, x2)

f <- function(x1, x2){cos(x1*x2)}

(2)

# The Second Taylor approximation of 
#f(x1, x2) = cos(x1, x2) is h(x1, x2) = 1-((pi^2)/8)*(x1^2)

h <- function(x1, x2){1-((pi^2)/8)*(x1^2)}

(3), (4)

library(rgl)
options(rgl.useNULL = TRUE)

x1 <- seq(-pi/4, pi/4, length.out=30)
x2 <- seq(pi/4, 3*pi/4, length.out=30)

my_f <- outer(x1, x2, FUN=f)
persp3d(x1, x2, my_f, col="red", xlab="x1", ylab="x2", shade=0.1, alpha=0.5,
main="Taylor approximation of cos(x1x2)")
#rglwidget(controllers = ) # This is to get it to print in pdf

my_h <- outer(x1, x2, FUN=h)
persp3d(x1, x2, my_h, col="blue", alpha = 0.5, add = TRUE)
rglwidget(controllers = ) # This is to get it to print in pdf

(5), (6)

error <- abs(my_f-my_h)

persp3d(x1, x2, error, col = "green", alpha = 0.5, xlab="x1", ylab="x2",
main="The error in second order Taylor expansion of cos(x1x2)")
rglwidget(controllers = ) # This is to get it to print in pdf

(7)

filled.contour(x1, x2, error)
grid(nx = 30, ny = 30)

Regarding the magnitude of the error, We get the actual approximation and the error is zero around the point (0, 1.7), but the error increases as it goes to every corner of the plot.

Problem 2

(a)

if we set \[\Sigma^{-1}=\begin{bmatrix} a & b \\ b & c \\ \end{bmatrix} \] then \[\Sigma^{-1}= \frac{1}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}}\begin{bmatrix} \sigma_{22} & -\sigma_{12} \\ -\sigma_{21} & \sigma_{11} \\ \end{bmatrix}=\begin{bmatrix} \frac{\sigma_{22}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}} & -\frac{\sigma_{12}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}} \\ -\frac{\sigma_{21}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}} & \frac{\sigma_{11}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}} \\ \end{bmatrix}=\begin{bmatrix} a & b \\ b & c \\ \end{bmatrix} \] so \[a= \frac{\sigma_{22}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}}\] \[b=-\frac{\sigma_{12}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}}\] \[c=\frac{\sigma_{11}}{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}}\] Also, \[\Sigma=\frac{1}{ac-b^2}\begin{bmatrix} c & -b \\ -b & a \\ \end{bmatrix}=\begin{bmatrix} c/ac-b^2 & -b/ac-b^2 \\ -b/ac-b^2 & a/ac-b^2 \\ \end{bmatrix}=\begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \\ \end{bmatrix}\] and \[|\Sigma|^{-1/2}=\sqrt{ac-b^2}=\frac{1}{\sqrt{\sigma_{11}\sigma_{22}-\sigma_{12}\sigma_{21}}}\]

Therefore, the function f(x) is \[f(x)=(2\pi)^{-1}|\Sigma|^{-1/2}exp\{-\frac{1}{2}ax_1^2+(a\mu_1+b\mu_2)x_1-\frac{1}{2}cx_2^2+(b\mu_1+c\mu_2)x_2-bx_1x_2-b\mu_1\mu_2-\frac{1}{2}a\mu_1^2-\frac{1}{2}c\mu_2^2\}\]

The Second order Taylor expansion for f(X) is \[f(x)=(2\pi)^{-1}|\Sigma|^{-1/2}\{1-ax_1^2+2(a\mu_1+b\mu_2)x_1-cx_2^2+2(b\mu_1+c\mu_2)x_2-2bx_1x_2-a\mu_1^2-2b\mu_1\mu_2-c\mu_2^2\}\]

(b), (c), (d) for case 1

a <- 1/0.91
b <- 0.3/0.91
c <- 1/0.91

f <- function(x1, x2){((2*(pi))^-1)*sqrt(a*c-(b^2))*exp((-1/2)*a*(x1^2)
                      +(-1/2)*c*(x2^2)-b*x1*x2)}
x1 <- seq(-1, 1, 0.1) 
x2 <- seq(-1, 1, 0.1)

my_f <- outer(x1, x2, FUN=f)
persp3d(x1, x2, my_f, col = "blue",shade = 0.1, alpha = 0.5, xlab="x1", ylab="x2")
#rglwidget(controllers = ) # This is to get it to print in pdf

g <- function(x1, x2){((2*(pi))^-1)*sqrt(a*c-(b^2))*(1-a*(x1^2)-c*(x2^2)
                      -2*b*x1*x2)}

my_g <- outer(x1, x2, FUN=g) 
persp3d(x1, x2, my_g, col = "red", add = TRUE, alpha = 0.5)
rglwidget(controllers = ) # This is to get it to print in pdf
error1 <- abs(my_f - my_g)
persp3d(x1, x2, error1, col = "green", alpha = 0.5)
rglwidget(controllers = ) # This is to get it to print in pdf
A <- matrix(c(1,-0.3, -0.3, 1), 2, 2)
ev_1 <- eigen(A)
ev_1
## eigen() decomposition
## $values
## [1] 1.3 0.7
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.7071068 -0.7071068
## [2,]  0.7071068 -0.7071068
filled.contour(x1, x2, my_f, plot.axes = { 
  axis(1)
  axis(2)
  lines(c(0, ev_1$vectors[1,1]*(ev_1$values[1])), 
        c(0, ev_1$vectors[2,1]*(ev_1$values[1])), col="red", lwd=2)
  lines(c(0, ev_1$vectors[1,2]*(ev_1$values[2])), 
        c(0, ev_1$vectors[2,2]*(ev_1$values[2])), col="blue", lwd=2)
})
grid(nx = 30, ny = 30)

(b), (c), (d) for case 2

a <- 1/0.36
b <- (-0.8)/0.36
c <- 1/0.36

f <- function(x1, x2){((2*(pi))^-1/sqrt(a*c-b*b))*exp((-1/2)*a*(x1^2)
                      +(-1/2)*c*(x2^2)-b*x1*x2)}
x1 <- seq(-1, 1, 0.1) 
x2 <- seq(-1, 1, 0.1)

my_f <- outer(x1, x2, FUN=f)
persp3d(x1, x2, my_f, col = "blue",shade = 0.1, alpha = 0.5, xlab="x1", ylab="x2")
#rglwidget(controllers = ) # This is to get it to print in pdf

g <- function(x1, x2){((2*(pi))^-1/sqrt(a*c-b*b))*(1-a*(x1^2)-c*(x2^2)
                      -2*b*x1*x2)}
my_g <- outer(x1, x2, FUN=g) 
persp3d(x1, x2, my_g, col = "red", add = TRUE, alpha = 0.5)
rglwidget(controllers = ) # This is to get it to print in pdf
error2 <- abs(my_f - my_g)
persp3d(x1, x2, error2, col = "green", alpha = 0.5, xlab="x1", ylab="x2")
rglwidget(controllers = ) # This is to get it to print in pdf
B <- matrix(c(1, 0.8, 0.8, 1), 2, 2)
ev_2 <- eigen(B)
ev_2
## eigen() decomposition
## $values
## [1] 1.8 0.2
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068
filled.contour(x1, x2, my_f, plot.axes = { 
  axis(1)
  axis(2)
  lines(c(0, ev_2$vectors[1,1]*(ev_2$values[1])), 
        c(0, ev_2$vectors[2,1]*(ev_2$values[1])), col="red", lwd=2)
  lines(c(0, ev_2$vectors[1,2]*(ev_2$values[2])), 
        c(0, ev_2$vectors[2,2]*(ev_2$values[2])), col="blue", lwd=2)
})
grid(nx = 30, ny = 30)

(c): Explanation

In both cases, the shapes of the constant value contours are \(ellipses\) and (0,0) is the center of value contours.

(d): Explanation

As can be seen from the eigenvectors superimposed in the contour plot, the eigenvectors are the major axes of the ellipses, and they indicate the direction of the principal components. The magnitude of these vectors equals the corresponding eigenvalue.