Expectation

# Expectation equals Mean when sample size approaches infinity.

# Set up a 3-row plotting area and adjust margins for better visibility.
par(mfrow = c(3, 1), mar = c(3, 4, 1, 1))
set.seed(12345678)
# Case 1: Small sample size (n = 10)
# The expectation of a chi-square distribution equals the degrees of freedom (df = 5).
x <- rchisq(n = 10, df = 5)  # Generate 10 random values from chi-square distribution
hist(x, 
     main = "Chi-Square Distribution (n=10, df=5)", 
     xlab = "Value", 
     col = "lightblue", 
     border = "black")

abline(v = mean(x), 
       col = "red", 
       lwd = 2)  # Add vertical line for sample mean

legend("topright", 
       legend = paste("Mean =", round(mean(x), 2)), 
       col = "red", 
       lwd = 2)

# Case 2: Larger sample size (n = 100)
x <- rchisq(n = 100, df = 5)  # Generate 100 random values
hist(x, 
     main = "Chi-Square Distribution (n=100, df=5)", 
     xlab = "Value", 
     col = "lightblue", 
     border = "black")

abline(v = mean(x), 
       col = "red", 
       lwd = 2)  # Add vertical line for sample mean

legend("topright", 
       legend = paste("Mean =", round(mean(x), 2)), 
       col = "red", 
       lwd = 2)

# Case 3: Even larger sample size (n = 10,000)
x <- rchisq(n = 10000, df = 5)  # Generate 10,000 random values
hist(x, 
     main = "Chi-Square Distribution (n=10,000, df=5)", 
     xlab = "Value", 
     col = "lightblue", 
     border = "black")

abline(v = mean(x), 
       col = "red", 
       lwd = 2)  # Add vertical line for sample mean

legend("topright", 
       legend = paste("Mean =", round(mean(x), 2)), 
       col = "red", 
       lwd = 2)

# Conclusion: As the sample size increases, the sample mean converges to the expected value (df = 5).

Variance

Variance is a measure of how spread out a set of values is. The formula for variance differs based on whether you are working with a population or a sample.

Population Variance (\(\sigma^2\))

The formula for population variance is:

\[ \sigma^2 = \frac{\sum_{i=1}^{N} (x_i - \mu)^2}{N} \]

Where:

  • \(x_i\) = each data point
  • \(\mu\) = population mean
  • \(N\) = total number of data points in the population

Sample Variance (\(s^2\))

The formula for sample variance is:

\[ s^2 = \frac{\sum_{i=1}^{n} (x_i - \bar{x})^2}{n - 1} \]

Where:

  • \(x_i\) = each data point
  • \(\bar{x}\) = sample mean
  • \(n\) = total number of data points in the sample
  • \(n - 1\) = degrees of freedom (used to correct for bias in the sample)

Example in R

Below is an example of how to compute sample variance in R:

# Sample data
data <- c(10, 12, 23, 23, 16, 19, 20, 22, 21, 18)

# Calculate population variance manually
population_variance <- sum((data - mean(data))^2) / length(data)
cat("Population Variance:", population_variance, "\n")
## Population Variance: 18.24
# Calculate sample variance manually
sample_variance <- sum((data - mean(data))^2) / (length(data)-1)
cat("Sample Variance:", sample_variance, "\n")
## Sample Variance: 20.26667
# Calculate sample variance automatically
sam_variance <- var(data)
cat("Sam Variance:", sam_variance, "\n")
## Sam Variance: 20.26667

Variance Calculation

set.seed(12345678)
# Set sample size
n <- 10  

# Generate 'n' random values from a normal distribution with mean = 100 and standard deviation = 5
x <- rnorm(n, mean = 100, sd = 5)  

# Plot histogram of the generated sample
# hist(x, 
#      main = "Histogram of Normally Distributed Data (n=10, mean=100, sd=5)", 
#      xlab = "Value", 
#      col = "lightblue", 
#      border = "black")

# Display the minimum and maximum values of the sample
range_values <- c(min(x), max(x))
cat("Minimum value:", range_values[1], "\n")
## Minimum value: 88.10468
cat("Maximum value:", range_values[2], "\n")
## Maximum value: 112.76
# Calculate variance using the formula: sum of squared differences divided by (n-1)
manual_variance <- sum((x - mean(x))^2) / (n - 1)
cat("Manual Variance Calculation:", manual_variance, "\n")
## Manual Variance Calculation: 51.56503
# Calculate standard deviation from the manually computed variance
manual_sd <- sqrt(manual_variance)
cat("Manual Standard Deviation Calculation:", manual_sd, "\n")
## Manual Standard Deviation Calculation: 7.180879
# Compare with built-in var() function in R
built_in_variance <- var(x)
cat("Built-in Variance Function Result:", built_in_variance, "\n")
## Built-in Variance Function Result: 51.56503
# Conclusion: The manually calculated variance should match the built-in var() function output.

Expectation and variance of linear function of random variables

Mathematical Properties of Random Variables

  1. Linear Transformation with a Scalar

If \(y = a \cdot x\), where \(a\) is a scalar, then:

  • The expected value of \(y\) is:

\[ E(y) = a \cdot E(x) \]

  • The variance of \(y\) is:

\[ Var(y) = a^2 \cdot Var(x) \]

  1. Additive Transformation

If \(y = x + a\), then:

  • The expected value of \(y\) is:

\[ E(y) = E(x) + a \]

  • The variance of \(y\) is:

\[ Var(y) = Var(x) \]

set.seed(12345678)
# Number of observations
n <- 10000

# Degrees of freedom for Chi-square distribution
df <- 10

# Generate n Chi-square distributed random variables with df degrees of freedom
x <- rchisq(n, df)

# Calculate and print the mean of x
mean(x)
## [1] 9.977369
# Calculate and print the variance of x
var(x)
## [1] 20.07793
# Create a new variable y which is 5 times x
y <- 5 * x

# Calculate and print the mean of y
mean(y)
## [1] 49.88684
# Calculate and print the variance of y
var(y)
## [1] 501.9481
# Calculate the mean of y using the linearity of expectation
5 * mean(x)
## [1] 49.88684
# Calculate the variance of y using the properties of variance
5^2 * var(x)
## [1] 501.9481
# Create a new variable z which is x plus 5
z <- 5 + x

# Calculate and print the mean of z
mean(z)
## [1] 14.97737
# Calculate and print the variance of z
var(z)
## [1] 20.07793
# Calculate the mean of z using the linearity of expectation
5 + mean(x)
## [1] 14.97737
# The variance of z remains the same as the variance of x
var(x)
## [1] 20.07793

Covariance Analysis

set.seed(12345678)
# Number of observations
n <- 10000

# Generate n Poisson distributed random variables with lambda = 100
x <- rpois(n, 100)

# Generate n Chi-square distributed random variables with df = 5
y <- rchisq(n, 5)

# Generate n t-distributed random variables with df = 100
z <- rt(n, 100)

# Set up plotting area to display three plots in one column
par(mfrow = c(3, 1), mar = c(3, 4, 1, 1))

# Plot x against y
plot(x, y)

# Plot x against z
plot(x, z)

# Plot y against z
plot(y, z)

# Calculate and print the variance of x
var(x)
## [1] 98.24031
# Calculate and print the variance of y
var(y)
## [1] 9.745051
# Calculate and print the variance of z
var(z)
## [1] 1.024255
# Calculate and print the covariance between x and y
cov(x, y)
## [1] -0.06837492
# Calculate and print the covariance between x and z
cov(x, z)
## [1] -0.2194389
# Calculate and print the covariance between y and z
cov(y, z)
## [1] 0.04445237
# Calculate covariance manually
sum((x-mean(x))*(y-mean(y)))/(n-1)
## [1] -0.06837492
sum((x-mean(x))*(z-mean(z)))/(n-1)
## [1] -0.2194389
sum((y-mean(y))*(z-mean(z)))/(n-1)
## [1] 0.04445237
set.seed(12345678)
# Define the number of observations
n <- 10000  

# Generate random variables with different distributions
a <- rnorm(n, mean = 100, sd = 5)  # Common base variable
x <- a + rpois(n, lambda = 100)    # Poisson-distributed component
y <- a + rchisq(n, df = 5)         # Chi-square distributed component
z <- a + rt(n, df = 100)           # t-distributed component

# Set up a 3-row plotting area and adjust margins for better visualization
par(mfrow = c(3, 1), mar = c(3, 4, 1, 1))

# Scatter plots to visualize relationships between variables
plot(x, y, 
     main = "Scatter Plot of X vs Y", 
     xlab = "X Values", 
     ylab = "Y Values", 
     col = "blue", 
     pch = 16)

plot(x, z, 
     main = "Scatter Plot of X vs Z", 
     xlab = "X Values", 
     ylab = "Z Values", 
     col = "red", 
     pch = 16)

plot(y, z, 
     main = "Scatter Plot of Y vs Z", 
     xlab = "Y Values", 
     ylab = "Z Values", 
     col = "green", 
     pch = 16)

# Calculate variances of individual variables
x_var <- var(x)
y_var <- var(y)
z_var <- var(z)

cat("Variance of X:", x_var, "\n")
## Variance of X: 123.4008
cat("Variance of Y:", y_var, "\n")
## Variance of Y: 35.17484
cat("Variance of Z:", z_var, "\n")
## Variance of Z: 25.7969
# Compute pairwise covariances between the variables
xy_cov <- cov(x, y)
xz_cov <- cov(x, z)
yz_cov <- cov(y, z)

cat("Covariance between X and Y:", xy_cov, "\n")
## Covariance between X and Y: 25.36729
cat("Covariance between X and Z:", xz_cov, "\n")
## Covariance between X and Z: 24.41266
cat("Covariance between Y and Z:", yz_cov, "\n")
## Covariance between Y and Z: 24.89724
# Combine variables into a matrix for further analysis
W <- cbind(x, y, z)

# Check the dimension of the matrix
cat("Dimension of W:", dim(W), "\n")
## Dimension of W: 10000 3
# Compute the covariance matrix of the combined variables
cov_matrix <- cov(W)
cat("Covariance Matrix:\n")
## Covariance Matrix:
print(cov_matrix)
##           x        y        z
## x 123.40081 25.36729 24.41266
## y  25.36729 35.17484 24.89724
## z  24.41266 24.89724 25.79690
# Compute the variance-covariance matrix (should be the same as covariance matrix)
var_matrix <- var(W)
cat("Variance-Covariance Matrix:\n")
## Variance-Covariance Matrix:
print(var_matrix)
##           x        y        z
## x 123.40081 25.36729 24.41266
## y  25.36729 35.17484 24.89724
## z  24.41266 24.89724 25.79690
# Compute the correlation matrix
cor_matrix <- cor(W)
cat("Correlation Matrix:\n")
## Correlation Matrix:
print(cor_matrix)
##           x         y         z
## x 1.0000000 0.3850339 0.4326856
## y 0.3850339 1.0000000 0.8265163
## z 0.4326856 0.8265163 1.0000000
# Conclusion: Covariance and correlation matrices help in understanding the relationships between variables.

Matrix Manipulations

# Create two matrices
a <- matrix(seq(10, 60, 10), nrow = 2, ncol = 3)
b <- matrix(seq(1, 6), nrow = 2, ncol = 3)

# Display matrices
cat("Matrix A:\n")
## Matrix A:
print(a)
##      [,1] [,2] [,3]
## [1,]   10   30   50
## [2,]   20   40   60
cat("Matrix B:\n")
## Matrix B:
print(b)
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
# Matrix addition
cat("Matrix Addition (A + B):\n")
## Matrix Addition (A + B):
print(a + b)
##      [,1] [,2] [,3]
## [1,]   11   33   55
## [2,]   22   44   66
# Matrix subtraction
cat("Matrix Subtraction (A - B):\n")
## Matrix Subtraction (A - B):
print(a - b)
##      [,1] [,2] [,3]
## [1,]    9   27   45
## [2,]   18   36   54
# Element-wise multiplication (dot product)
cat("Element-wise Multiplication (A * B):\n")
## Element-wise Multiplication (A * B):
print(a * b)
##      [,1] [,2] [,3]
## [1,]   10   90  250
## [2,]   40  160  360
# Element-wise division
cat("Element-wise Division (A / B):\n")
## Element-wise Division (A / B):
print(a / b)
##      [,1] [,2] [,3]
## [1,]   10   10   10
## [2,]   10   10   10

Matrix multiplication

c <- matrix(c(1, 1, 1, 1, 3, 4, 2, 1, 30, 50, 40, 25), nrow = 4, ncol = 3)
d <- matrix(c(20000, 10000, 1000, 1000, 300, 20), nrow = 3, ncol = 2)

# Perform matrix multiplication
t <- c %*% d
cat("Matrix Multiplication Result (C %*% D):\n")
## Matrix Multiplication Result (C %*% D):
print(t)
##        [,1] [,2]
## [1,]  80000 2500
## [2,] 110000 3200
## [3,]  80000 2400
## [4,]  55000 1800

Matrix Inverse

# Matrix inverse (only for square matrices)
t <- t[1:2, 1:2]  # Extract a 2x2 sub-matrix
ti <- solve(t)    # Compute the inverse

cat("Inverse of Matrix T:\n")
## Inverse of Matrix T:
print(ti)
##               [,1]          [,2]
## [1,] -0.0001684211  0.0001315789
## [2,]  0.0057894737 -0.0042105263
# Verify inverse property: T * T⁻¹ = I
cat("Verification (T⁻¹ * T):\n")
## Verification (T⁻¹ * T):
print(ti %*% t)
##              [,1]         [,2]
## [1,] 1.000000e+00 5.551115e-17
## [2,] 5.684342e-14 1.000000e+00
cat("Verification (T * T⁻¹):\n")
## Verification (T * T⁻¹):
print(t %*% ti)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

Matrix transpose

c <- matrix(c(1, 1, 1, 4, 30, 50), nrow = 2, ncol = 3)
cat("Matrix C:\n")
## Matrix C:
print(c)
##      [,1] [,2] [,3]
## [1,]    1    1   30
## [2,]    1    4   50
cat("Transpose of Matrix C:\n")
## Transpose of Matrix C:
print(t(c))
##      [,1] [,2]
## [1,]    1    1
## [2,]    1    4
## [3,]   30   50
# Properties of transpose
A <- matrix(c(1, 1, 1, 4, 30, 50), nrow = 2, ncol = 3)
B <- matrix(c(1000, 300, 20, 20000, 10000, 1000), nrow = 3, ncol = 2)

# Transpose properties verification
cat("Transpose of (A %*% B):\n")
## Transpose of (A %*% B):
print(t(A %*% B))
##       [,1]   [,2]
## [1,]  1900   3200
## [2,] 60000 110000
cat("Transpose of B %*% Transpose of A:\n")
## Transpose of B %*% Transpose of A:
print(t(B) %*% t(A))
##       [,1]   [,2]
## [1,]  1900   3200
## [2,] 60000 110000

Special Matrices

  • Symmetric Matrix: A matrix \(A\) is symmetric if it is equal to its transpose, i.e., \(A = A^T\).
  • Diagonal Matrix: A matrix in which all off-diagonal elements are zero, meaning only the diagonal elements can be nonzero.
  • Identity Matrix: A special diagonal matrix where all diagonal elements are 1 and all off-diagonal elements are 0.
  • Orthogonal Matrix: A matrix \(A\) is orthogonal if the product of the matrix and its transpose equals the identity matrix, i.e., \(A A^T = I\).
  • Singular Matrix: A square matrix that does not have an inverse, meaning its determinant is zero.
# Define matrices
A <- matrix(c(1, 2, 2, 1), 2, 2)  # Symmetric matrix
D <- diag(c(3, 5, 7))  # Diagonal matrix
I <- diag(3)  # Identity matrix
Q <- matrix(c(0, 1, -1, 0), 2, 2)  # Orthogonal matrix
S <- matrix(c(2, 4, 1, 2), 2, 2)  # Singular matrix

A
##      [,1] [,2]
## [1,]    1    2
## [2,]    2    1
D
##      [,1] [,2] [,3]
## [1,]    3    0    0
## [2,]    0    5    0
## [3,]    0    0    7
I
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
Q
##      [,1] [,2]
## [1,]    0   -1
## [2,]    1    0
S
##      [,1] [,2]
## [1,]    2    1
## [2,]    4    2
# Check symmetry
is_symmetric <- all(A == t(A))
print(paste("Is A symmetric? ", is_symmetric))
## [1] "Is A symmetric?  TRUE"
# Check identity matrix
is_identity <- all(I %*% I == I)
print(paste("Is I an identity matrix? ", is_identity))
## [1] "Is I an identity matrix?  TRUE"
# Check orthogonality
is_orthogonal <- all(Q %*% t(Q) == diag(2)) # diag(2) creates a 2x2 identity matrix
print(paste("Is Q orthogonal? ", is_orthogonal))
## [1] "Is Q orthogonal?  TRUE"
# Check if matrix is singular (determinant should be zero)
is_singular <- det(S) == 0
print(paste("Is S singular? ", is_singular))
## [1] "Is S singular?  TRUE"

Rank of a Matrix

Rank of a matrix refers to the number of linearly independent rows or columns in the matrix. It indicates the dimension of the column space or row space, which represents the number of independent directions or basis vectors in the matrix.

the qr() function performs QR decomposition (also called QR factorization) of a matrix. QR decomposition is a technique used in numerical linear algebra to factorize a given matrix into an orthogonal (or unitary) matrix and an upper triangular matrix, which can be useful for solving linear systems, computing matrix rank, and performing least squares regression.

# Define a matrix A
A <- matrix(c(1, 2, 3, 
              4, 5, 6, 
              7, 8, 9), nrow = 3, byrow = TRUE)

# Display the matrix
print(A)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9
# Compute the rank of matrix A
rank_A <- qr(A)$rank
print(paste("Rank of matrix A:", rank_A))
## [1] "Rank of matrix A: 2"
# The matrix has dependent rows, so its rank is 2.

Trace of a Matrix

Trace of a square matrix is the sum of its diagonal elements. It represents the sum of the eigenvalues of the matrix and has important applications in linear algebra, such as measuring matrix similarity and transformations.

# Define a square matrix B
B <- matrix(c(4, 2, 1, 
              0, 5, 3, 
              7, 8, 6), nrow = 3, byrow = TRUE)

# Display the matrix
print(B)
##      [,1] [,2] [,3]
## [1,]    4    2    1
## [2,]    0    5    3
## [3,]    7    8    6
# Compute the trace of the matrix (sum of diagonal elements)
trace_B <- sum(diag(B))
print(paste("Trace of matrix B:", trace_B))
## [1] "Trace of matrix B: 15"