Probability and Statistics (UMA401)
Experiment 6: Joint probability mass and density functions
# Load the required package
library(pracma)
(1) The joint probability density of two random variables X and Y is
2(2x + 3y)/5; 0 ≤ x, y ≤ 1
F(x) =
0; elsewhere
Write a R-code to
(i) check that it is a joint density function or not? (Use integral2()).
f <- function(x, y) 2*(2*x + 3*y)/5
total_prob <- integral2(f, 0, 1, 0, 1)$Q
is_joint_density <- abs(total_prob - 1) < 1e-6
(ii) find marginal distribution g(x) at x = 1.
g_x <- function(y) f(1, y)
marginal_g <- integral(g_x, 0, 1)
(iii) find the marginal distribution h(y) at y = 0.
h_y <- function(x) f(x, 0)
marginal_h <- integral(h_y, 0, 1)
(iv) find the expected value of g(x, y) = xy.
expected_value <- integral2(function(x,y) x*y*f(x,y), 0, 1, 0, 1)$Q
(2) The joint probability mass function of two random variables X and Y is
f (x, y) = (x + y)/30; x = 0, 1, 2, 3; y = 0, 1, 2
Write a R-code to
# Define possible values
x_vals <- 0:3
y_vals <- 0:2
(i) display the joint mass function in rectangular (matrix) form.
joint_pmf <- matrix(0, nrow = length(x_vals), ncol = length(y_vals))
for (i in 1:length(x_vals)) {
for (j in 1:length(y_vals)) {
x <- x_vals[i]
1
Probability and Statistics (UMA401)
y <- y_vals[j]
joint_pmf[i, j] <- (x + y) / 30
cat("Joint PMF Matrix:\n")
print(joint_pmf)
OR
# Create joint PMF matrix
joint_pmf <- outer(x_vals, y_vals, function(x, y) (x + y)/30)
# Optional: name rows and columns for readability
rownames(joint_pmf) <- paste0("x=", x_vals)
colnames(joint_pmf) <- paste0("y=", y_vals)
print(joint_pmf)
(ii) check that it is joint mass function or not? (use: Sum())
sum_pmf <- sum(joint_pmf)
cat("Sum of PMF =", sum_pmf, "\n") # should be 1
(iii) find the marginal distribution g(x) for x = 0, 1, 2, 3. (Use:apply())
g_x <- apply(joint_pmf, 1, sum)
cat("Marginal g(x):\n")
print(g_x)
(iv) find the marginal distribution h(y) for y = 0, 1, 2. (Use:apply())
h_y <- apply(joint_pmf, 2, sum)
cat("Marginal h(y):\n")
print(h_y)
2
Probability and Statistics (UMA401)
(v) find the conditional probability at x = 0 given y = 1.
# # P(X=0, Y=1) = joint_pmf[1, 2] (X=0 is row 1, Y=1 is column 2)
numerator <- joint_pmf[1, 2]
# P(Y=1) = sum of column 2
denominator <- sum(joint_pmf[, 2])
P_x0_given_y1 <- numerator / denominator
cat("P(X = 0 | Y = 1) =", P_x0_given_y1, "\n")
(vi) find E(x), E(y), E(xy), Var(x), Var(y), Cov(x, y) and its correlation coefficient.
# pxy <- [Link](joint_pmf)
# Create x and y grid
x_grid <- rep(x_vals, each = length(y_vals)) # 0 0 0 1 1 1 ...
y_grid <- rep(y_vals, times = length(x_vals)) # 0 1 2 0 1 2 ...
# E[X], E[Y], E[XY]
E_x <- sum(x_grid * pxy)
E_y <- sum(y_grid * pxy)
E_xy <- sum(x_grid * y_grid * pxy)
# Var(X), Var(Y)
Var_x <- sum((x_grid - E_x)^2 * pxy)
Var_y <- sum((y_grid - E_y)^2 * pxy)
# Cov(X, Y)
Cov_xy <- E_xy - E_x * E_y
3
Probability and Statistics (UMA401)
# Correlation coefficient
Corr_xy <- Cov_xy / sqrt(Var_x * Var_y)
# Print results
cat("E(X) =", E_x, "\n")
cat("E(Y) =", E_y, "\n")
cat("E(XY) =", E_xy, "\n")
cat("Var(X) =", Var_x, "\n")
cat("Var(Y) =", Var_y, "\n")
cat("Cov(X,Y) =", Cov_xy, "\n")
cat("Correlation Coefficient =", Corr_xy, "\n")
4
5