Individual assignment III

Task 1

From the theoretical examples, I decided to solve the Exercise 2: \ Consider the random vector \(\mathbf{X}\) with the two-dimensional normal distribution \(\mathcal{N}_2\left(\begin{pmatrix} 1 \\ 2 \end{pmatrix}, \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\right)\) and the conditional random vector \(\mathbf{Y}|\mathbf{X}\sim \mathcal{N}_2 \left(\begin{pmatrix} X_1 \\ X_1 + X_2 \end{pmatrix}, \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\right)\).

Determine the distribution of \(Y_2|Y_1\) and \(\mathbf{W} = \mathbf{X}- \mathbf{Y}.\)

The first step is going to be to find the distribution of \((\mathbf{X}^T, \mathbf{Y}^T)^T\). From the Exercise 1, we obtain that \[\begin{align*} (\mathbf{X}^T, \mathbf{Y}^T)^T \sim \mathcal{N}_4(\boldsymbol{\mu}, \Sigma), \end{align*}\] where \(\boldsymbol{\mu} = (1,2,1,3)^T\) and \(\Sigma = \begin{pmatrix} 2 & 1 & 2 & 3 \\ 1 & 2 & 1 & 3 \\ 2 & 1 & 3 & 3 \\ 3 & 3 & 3 & 7 \end{pmatrix}\)

To find the conditional distribution \(Y_2|Y_1\), we easily derive the distribution of \((Y_1, Y_2)^T:\) It follows from the joint distribution of \((\mathbf{X}, \mathbf{Y})^T\) that \(\mathbf{Y}\sim \mathcal{N}_2 \left(\begin{pmatrix} 1 \\ 3 \end{pmatrix}, \begin{pmatrix} 3 & 3 \\ 3 & 7 \end{pmatrix}\right)\). Thus, according to the computations on the lab session, \(Y_2|Y_1 \sim \mathcal{N}(3 + (Y_1 - 1)\frac{3}{3}, 7 - \frac{9}{3}) = \mathcal{N}(2 + Y_1, 4).\)

Regarding the distribution of \(\mathbf{W}\), one should realise that a linear transformation preserves normality. From \((\mathbf{X}, \mathbf{Y})^T\) we get \(\mathbf{X}- \mathbf{Y}\) by a linear transformation given by a matrix \(A = \begin{pmatrix} 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \end{pmatrix}\).

The expectation of \(\mathbf{X}- \mathbf{Y}\) is computed easily as \((0,-1)^T\) and the variance matrix is given as follows: \(\mathsf{var}(\mathbf{X}- \mathbf{Y}) = A ~ \mathsf{var}(\mathbf{X}^T, \mathbf{Y}^T) A^T\). Thus, \(\mathbf{W}\sim \mathcal{N}_2 \left(\begin{pmatrix} 0 \\ -1 \end{pmatrix}, \begin{pmatrix} 1 & 0 \\ 0 & 3 \end{pmatrix}\right).\)  

Now, to look at the distribution of \(\mathbf{W}= \mathbf{X}- \mathbf{Y}\), we generate two random samples:

library("MASS")
##two random samples from the distributions we consider
set.seed(123)
X.Y = mvrnorm(1000, mu = c(1,2,1,3), Sigma = matrix(c(2,1,2,3,1,2,1,3,2,1,3,3,3,3,3,7),4))


W = X.Y[, c(1,2)] - X.Y[, c(3,4)]
mean(W[,1]) ##should be around 0 

## [1] 0.006547564

mean(W[,2]) ##should be around 1

## [1] -0.9550145

## two dimensional kernel density estimate
bivn.kde1 <- kde2d(W[,1], W[,2], n = 50)

par(mfrow = c(1,2))
persp(bivn.kde1, phi = 45, theta = 30, shade = .1, border = NA)
##now compare it with the expected distribution
N = mvrnorm(1000, mu = c(0,-1), Sigma = matrix(c(1,0,0,3), 2))
bivn.kde2 <- kde2d(N[,1], N[,2], n = 50)
persp(bivn.kde2, phi = 45, theta = 30, shade = .1, border = NA)

Next, we can plot a scatterplot of \(\mathbf{W}\) to see if the data corresponds with parameters of the theoretical disribution of \(\mathbf{W}\) which we derived:

plot(W[,2] ~ W[,1], pch = 21, col = "gray")
abline(lm(W[,2] ~ W[,1]), col  = "red", lwd = 2)

plot(N[,2] ~ N[,1], pch = 21, col = "gray")
abline(lm(W[,2] ~ W[,1]), pch = 21, col = "red")

Task 2

Let \(X \sim \mathcal{N}(0, 1), ~ W \sim Alt(1/2)\) be independent. Let us define a random variable \[\begin{equation}\label{Y} Y = X \mathbf{1}_{[W = 0]} - X \mathbf{1}_{[W = 1]} = X(\mathbf{1}_{[W = 0]} - \mathbf{1}_{[W = 1]}) \end{equation}\]

We can compute the distribution function of \(Y\): \[\begin{align*} F_Y (y) &= \mathsf{P}[X(\mathbf{1}_{[W=0]} - \mathbf{1}_{[W=1]}) \leq y] \\ &= \mathsf{P}[X(\mathbf{1}_{[W=0]} - \mathbf{1}_{[W=1]}) \leq y ~ | ~ W=0] \cdot \mathsf{P}[W=0] + \mathsf{P}[X(\mathbf{1}_{[W=0]} - \mathbf{1}_{[W=1]}) \leq y ~ | ~ W=1] \cdot \mathsf{P}[W=1] \\ &= \mathsf{P}[X \leq y ~ | ~ W = 0] \cdot \frac{1}{2} + \mathsf{P}[-X \leq y ~ | ~ W=1] \cdot \frac{1}{2} \\ &= \frac{1}{2} \cdot \mathsf{P}[X \leq y] + \frac{1}{2} \cdot \mathsf{P}[-X \leq y] \\ &= \Phi(y) \end{align*}\]

Thus, \(Y \sim \mathcal{N}(0,1)\) as well.

Now, let us define a random variable \[\begin{align*} Z = X + Y = \begin{cases} 2X & \text{with prob.} ~ \frac{1}{2} \\ 0 & \text{with prob.} ~ \frac{1}{2} \end{cases} \end{align*}\] It means that \(Z\) is a mixed random variable with a density \[\begin{align*} f_Z(z) = \frac{1}{2} \delta_0(z) + \frac{1}{2}\frac{1}{\sqrt{8\pi}} \exp\left\{-\frac{z^2}{8}\right\} = \frac{1}{2} \delta_0(z) + \frac{1}{4\sqrt{2\pi}} \exp\left\{-\frac{z^2}{8}\right\} \end{align*}\]

Hence, \(Z\) is not normally distributed and it is a linear combination of \(X\) and \(Y\), so the vector \((X,Y)^T\) is not normally distributed.

From the following histograms, we can see that the random sample \(Y_1, \ldots, Y_{1000}\) of random variables defined as \(Y\) really comes from the standard normal distribution.

set.seed(123)
X = rnorm(1000, mean = 0, sd = 1)

W = rbinom(1000, 1, 0.5)

Y = NULL
for (i in 1:1000){if (W[i] == 0){Y = c(Y, X[i])}
                  else {Y = c(Y, -X[i])}}

par(mfrow = c(1,2))
hist(X, prob = TRUE)
curve(dnorm(x, mean = 0, sd = 1), add = TRUE)

hist(Y, prob = TRUE)
curve(dnorm(x, mean = 0, sd = 1), add=TRUE)

We can use histogram as well to show that \(X+Y\) is not normally distributed:

hist(X+Y, pch = 21, bg = "lightblue", prob = TRUE)
curve(dnorm(x, mean = mean(X+Y), sd = sd(X+Y)), add = TRUE)

Obviously, there is no normality there.

We can use another graphical tool to visualise that the normality of \((X,Y)^T\) is not fulfilled. First, we can plot the random sample \((X_1, Y_1)^T, \ldots, (X_{1000}, Y_{1000})^T:\)

plot(Y ~ X, pch = 21, bg = "lightblue")

We will add the joint normal density with corresponding parameters:

mu <- c(mean(X), mean(Y))   ## sample mean vector
sigma <- matrix(c(var(X), cov(X,Y), cov(X,Y), var(Y)), nrow = 2, byrow = TRUE)              

sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

ellipse <- function(s,t){
                    u<-c(s,t)-mu; 
                    e <- (-1) * (u %*% sigma.inv %*% u) / 2; 
                    exp(e)/(2 * pi * sqrt(det(sigma.inv)))}

n <- 60
x <- (0:(n-1))*2 - 50
y <- (0:(n-1))*2 - 50

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))

library(RColorBrewer)
plot(Y ~ X, pch = 21, bg = "lightblue")
contour(x,y,matrix(z,n,n), levels=(0:15), col = terrain.colors(16), add=TRUE, lwd = 1)