##Libraries

library(plot3D)

Problem 1a(i, ii, iii, iv)

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

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

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

f.out <- outer(x1vals, x2vals, f)

h.out <- outer(x1vals, x2vals, h)
persp3D(x1vals, x2vals, f.out, xlab = "x1", ylab = "x2", main = "Taylor approximation of cos(x1x2)", col = "red")

persp3D(x1vals, x2vals, h.out, add=TRUE, col = "blue")

Problem #1a(v, vi)

error <- function(x1, x2) {
    z = abs(f(x1,x2) - h(x1, x2))
}

error.out <- outer(x1vals, x2vals, error)
persp3D(x1vals, x2vals, error.out, xlab = "x1", ylab = "x2", main = "The error in the Taylor approximation of cos(x1x2)")

Problem 1a(vii)

contour(x1vals, x2vals, error.out, xlab = "x1", ylab = "x2", main="Contour plot of error of Taylor approximation of cos(x1x2)")

Problem 2a. Below we start with the bivariate density and we will compute partial derivatives

\[\begin{equation} f(x) = \frac{1}{2\pi\sqrt{det(\Sigma)}}*exp[-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)] \\ \textrm{Since } \Sigma \textrm{ is symmetric 2x2 matrix, } det(\Sigma) = \sigma_{11}\sigma_{22} - \sigma_{12}^2 \\ \textrm{ Let } c = \sigma_{11}\sigma_{22} - \sigma_{12}^2\\ \textrm{Due to unfamiliarity with Latex commands, much of the algebra will be omitted in the calculation of the partial derivatives.}\\ \textrm{If the reader wishes to see algebra, please contact the author at alzate.hunter@csu.fullerton.edu}\\ f_{x_1} = -\frac{f}{c}(\sigma_{22}(x_1 - \mu_1) - \sigma_{12}(x_2 - \mu_2))\\ f_{x_2} = -\frac{f}{c}(-\sigma_{12}(x_1 - \mu_1) + \sigma_{11}(x_2 - \mu_2))\\ f_{x_1x_1} = -\frac{f}{c}\sigma_{22} - \frac{1}{c}(\sigma_{22}(x_1 - \mu_1) - \sigma_{12}(x_2 - \mu_2))f_{x_1}\\ f_{x_2x_2} = -\frac{f}{c}\sigma_{11} - \frac{1}{c}(\sigma_{11}(x_2 - \mu_2) - \sigma_{12}(x_1 - \mu_1))f_{x_2}\\ f_{x_1x_2} = \frac{f}{c}\sigma_{12} - \frac{1}{c}(\sigma_{11}(x_2 - \mu_2) - \sigma_{12}(x_1 - \mu_1))f_{x_1}\\ \textrm{The 2nd order Taylor expansion of a function of 2 variables was derived in class. In that formula we evaluate the partial derivatives at the vector } \mathbf{x_0}.\\ \textrm{Conveniently the first partial deriviatives evaluate to 0 when we plug in this vector } \mathbf{x_0}.\\ \end{equation}\]

Simplifying the 2nd order expansion yields

\[\begin{eqnarray} f(\mathbf{x}) &\approx\ & f(\mu_1, \mu_2) + 0 + 0 + \frac{1}{2}f_{x_1x_1}(\mu_1,\mu_2)(x_1 - \mu_1)^2 +\frac{1}{2}f_{x_2x_2}(\mu_1,\mu_2)(x_2 - \mu_2)^2 + f_{x_1x_2}(\mu_1, \mu_2)(x_1 - \mu_1)(x_2-\mu_2)\\ & = & \frac{1}{2\pi\sqrt{c}} + \frac{-\sigma_{22}(x_1-\mu_1)^2}{4\pi*c^\frac{3}{2}} + \frac{-\sigma_{11}(x_2-\mu_2)^2}{4\pi*c^\frac{3}{2}}+ \frac{\sigma_{12}(x_1-\mu_1)(x_2-\mu_2)}{2\pi*c^\frac{3}{2}}\\ &=& \frac{1}{2\pi*\sqrt{c}}(1 - \frac{\sigma_{22}(x_1-\mu_1)^2}{2c} - \frac{\sigma_{11}(x_2-\mu_2)^2}{2c} + \frac{\sigma_{12}}{c}(x_1-\mu_1)(x_2-\mu_2))\\ \end{eqnarray}\]

Problem 2a in R code. The exact density and the 2nd order approximation.

density_f <- function(x1, x2) {
    expression <- (x1-mu[1])*(sigma_22*(x1-mu[1]) - sigma_12*(x2-mu[2])) + (x2-mu[2])*(-1*sigma_12*(x1-mu[1]) + sigma_11*(x2-mu[2]))
    z = ((2*pi)^(-1*p/2))*(1/det(S))*exp((-0.5*(1/c)) * expression)
    return(z)
}

approx_f <- function(x1, x2) {
    a = (1/(2*pi*(sqrt(c))))
    expansion = 1 - (sigma_22/(2*c))*((x1 - mu[1])^2) - (sigma_11/(2*c))*((x2 - mu[2])^2) + (sigma_12/c)*(x1-mu[1])*(x2-mu[2])
    z = a*expansion
    return(z)
}

Problem 2b(Case i.)

sigma_11 <- 1
sigma_12 <- -1*0.3
sigma_22 <- 1

S <- matrix(c(sigma_11, sigma_12, sigma_12, sigma_22), nrow=2, byrow=2)
mu <- c(0,0)
c <- det(S)
p <- 2

x1vals <- seq(-3, 3, length.out = 30)
x2vals <- seq(-3, 3, length.out = 30)

densityf_out <- outer(x1vals, x2vals, density_f)
approxf_out <- outer(x1vals, x2vals, approx_f)

persp3D(x1vals, x2vals, densityf_out, xlab = "x1", ylab = "x2", main = "Taylor approximation of Bivarite Normal Part (i)", col ="red")
persp3D(x1vals, x2vals, approxf_out, add=TRUE, col = "blue")

Problem 2cd(Case i)

error_Two <- function(x1, x2) {
    
    z = abs(density_f(x1,x2) - approx_f(x1, x2))
}
                   
errorTwo.out <- outer(x1vals, x2vals, error_Two)
persp3D(x1vals, x2vals, errorTwo.out, xlab = "x1", ylab = "x2", main = "The error in the second order Taylor expansion of Bivariate Normal with covariance = -0.3")

contour(x1vals, x2vals, errorTwo.out, xlab = "x1", ylab = "x2", main="Contour plot of error of Bivarite Normal with covariance = -0.3")

e <- eigen(S)
segments(0,0, e$vectors[1,1], e$vectors[2,1])
segments(0,0, e$vectors[1,2], e$vectors[2,2])

The shape of the constant contour values when the covariances = -0.3 forms ellipses. It appears the center of the constant contour values is the point (0,O)

The eigenvectors appears to be the relative angle formed by the major and minor axes of the ellipse formed by the constant contours. We will see if this is the case when the covariance = 0.8 in problem 2abcd(Case ii)

Problem 2abcd(Case ii)

sigma_11 <- 1
sigma_12 <- 0.8
sigma_22 <- 1
S <- matrix(c(sigma_11, sigma_12, sigma_12, sigma_22), nrow=2, byrow=2)
mu <- c(0,0)
c <- det(S)
p <- 2

densityf_out <- outer(x1vals, x2vals, density_f)
approxf_out <- outer(x1vals, x2vals, approx_f)

persp3D(x1vals, x2vals, densityf_out, xlab = "x1", ylab = "x2", main = "Taylor approximation of Bivariate Normal Part (Case ii)", col ="red")
persp3D(x1vals, x2vals, approxf_out, add=TRUE, col = "blue")

error_Two <- function(x1, x2) {
    
    z = abs(density_f(x1,x2) - approx_f(x1, x2))
}
                   
errorTwo.out <- outer(x1vals, x2vals, error_Two)
persp3D(x1vals, x2vals, errorTwo.out, xlab = "x1", ylab = "x2", main = "The error in the second order Taylor expansion of Bivariate Normal with covariance = 0.8")

contour(x1vals, x2vals, errorTwo.out, xlab = "x1", ylab = "x2", main="Contour plot of error of Bivarite Normal with covariance = 0.8")
e <- eigen(S)
segments(0,0, e$vectors[1,1], e$vectors[2,1])
segments(0,0, e$vectors[1,2], e$vectors[2,2])

It appears that again the eigenvectors form the major and minor axes of the ellipse formed by the contour values.