Final Year Solve 24-23

2024

1.(a) Create a vector x = (1,2,2,3,3,3,4,4,4,4,5,5,5,5,5) by using rep() and find multiplication of x.

x <- rep(1:5, c(1,2,3,4,5))
x

 [1] 1 2 2 3 3 3 4 4 4 4 5 5 5 5 5

1.(b) Create vector y = N * e^x / x!, where N = 1:length(x). Find sum of y.

N <- 1:length(x)
y <- N*exp(x)/factorial(x)
y

 [1]  2.718282  7.389056 11.083584 13.390358 16.737947 20.085537
 [7] 15.924460 18.199383 20.474306 22.749229 13.604540 14.841316
[13] 16.078092 17.314869 18.551645

y_sum <- sum(y)
y_sum

[1] 229.1426

1.(c) Find index positions where 10 < y < 15. Find the sum of values.

indices <- which(y>10 & y<15)
indices

[1]  3  4 11 12

val_sum <- sum(y[indices])
val_sum

[1] 52.9198

1.(d) Create int.y as nearest integer of y, count divisible by 4

int.y <- round(y)
int.y

 [1]  3  7 11 13 17 20 16 18 20 23 14 15 16 17 19

n <- length(which(int.y%%4==0))
n

[1] 4

1.(e) Function to create y and calculate mean

calc_mean <- function(x){
  N <- 1:length(x)
  y <- N*exp(x)/factorial(x)
  am <- mean(y)
  if (any(y<=0)){
    gm <- NA
    hm <- NA
  } else {
    gm <- exp(mean(log(y)))
    hm <- 1/mean(1/y)
  }
  return(list(AM=am,GM=gm,HM=hm))
}

1.(f) Evaluate the function for X.

means_result <- calc_mean(x)
cat("Arithmetic mean :", means_result$AM, "\n")

Arithmetic mean : 15.27617 

cat("Geometric mean  :", means_result$GM, "\n")

Geometric mean  : 13.90898 

cat("Harmonic mean   :", means_result$HM, "\n")

Harmonic mean   : 11.5125 

2.(a) Using seed 100, create three stdnormal random vectors of length 100 (x, x1,x2).

1. Create vectors y,z [y=x+x1, z=y+x2] and Evaluate sum((x-xbar)(y-ybar)), sum((x-xbar)(z-zbar)), sum((y-ybar)(z-zbar))

2. create matrix A(100x3) where each column is filled by each vector x,y,z. B=(A^T A) then B^-1 = ? Find sum of diagonal elements of B.

# (i)
set.seed(100)
x <- rnorm(100)
x1 <- rnorm(100)
x2 <- rnorm(100)

y <- x+x1
z <- y+x2
spxy <- sum((x-mean(x))*(y-mean(y)))
spxz <- sum((x-mean(x))*(z-mean(z)))
spxy

[1] 92.27519

spxz

[1] 103.4904

# (ii)
A <- matrix(c(x,y,z), byrow = FALSE, ncol = 3)
head(A)

            [,1]        [,2]       [,3]
[1,] -0.50219235 -0.83511570 -0.8069439
[2,]  0.13153117  1.49464487  1.1379415
[3,] -0.07891709 -0.54806443  0.3045619
[4,]  0.88678481  1.72966044  2.2430257
[5,]  0.11697127 -1.34102245 -0.3228195
[6,]  0.31863009 -0.08167583 -1.1031549

B <- t(A) %*% A
B

          [,1]      [,2]     [,3]
[1,] 103.14396  92.27928 103.4982
[2,]  92.27928 144.18476 133.6212
[3,] 103.49823 133.62118 230.0049

n <- nrow(B)
diag_sum <- 0
# sum(diag(B))
for (i in 1:n){
    diag_sum <- diag_sum + B[i,i]
}
diag_sum

[1] 477.3336

2.(b) Function to count entries > k in each column of a matrix.

count_gt <- function(mat, k) {
  apply(mat, 2, function(col) sum(col > k))
}

2.(c) Evaluate the func with (A,2)

count_gt(A,2)

[1] 3 4 8

3.(a) Gamma distribution calculations (shape=3, scale=0.25) (i) P(X ≥ 50) (ii) P(20 < X < 40) (iii) Find k such that P(X > k) = 0.1

#(i)
p50_ <- 1-pgamma(50, shape=3, scale=0.25)
p50_

[1] 0

#(ii)
p20_40 <- pgamma(40, shape=3, scale=0.25)-pgamma(20, shape=3, scale=0.25)
p20_40

[1] 0

#(iii)
k <- qgamma(0.9, shape = 3, scale = 0.25)  # Since P(X > k) = 0.1 => P(X ≤ k) = 0.9
k

[1] 1.33058

3.(b) Generate 5000 observations from gamma(3,0.25) [seed 101] (i) Estimate P(X ≥ 50) (ii) Estimate P(20 < X < 40) (iii) Find deciles (10%, 20%, …, 90%)

set.seed(101)
x<-rgamma(5000, shape=3, scale=0.25)
#(i)
est_px50_ <- mean(x>=50)
est_px50_

[1] 0

#(ii)
est_px20_50 <- mean(x>20 & x<40)
est_px20_50

[1] 0

#(iii)
deciles <- quantile(x, probs = seq(0.1, 0.9, by = 0.1))
deciles

      10%       20%       30%       40%       50%       60%       70% 
0.2838650 0.3900256 0.4802828 0.5736065 0.6662751 0.7708639 0.9059050 
      80%       90% 
1.0696102 1.3150904 

---

Question 4

(a) Finding the minimum value of the function [ l(x) = x^4 + 4x^3 + 2x^2 + 12x + 5] requires solving the transcendental equation [ f(x) = 0] where [ f(x) = 4x^3 + 12x^2 + 4x + 12] and [ f’(x) = 12x^2 + 24x + 4]

(i) Draw the functions ( l(x) ) and ( f(x) ) on the same graph for the interval [ -4 < x < 1]

(ii) Write a function to find the solution of the transcendental equation ( f(x) = 0 ). The function should return:

• the estimated solution,
• the value of the function at the estimated solution, and
• the number of iterations required.

---

(b) Write a function that sums the terms of the given series.

---

---

Workflow for Question (a)

Step 1: Define the functions

• Define the original function [ l(x) = x^4 + 4x^3 + 2x^2 + 12x + 5]
• Define its first derivative [ f(x) = l’(x) = 4x^3 + 12x^2 + 4x + 12]
• Define the derivative of ( f(x) ) [ f’(x) = 12x^2 + 24x + 4]

---

Step 2: Visualization

• Choose the interval ( -4 < x < 1 )

• Plot:

  • ( l(x) )
  • ( f(x) )

• Draw both curves on the same graph

• Use the plot to understand where the minimum of ( l(x) ) may occur

---

Step 3: Set up the optimization condition

• Recall that the minimum of ( l(x) ) occurs when: [ f(x) = 0]
• This equation cannot be solved easily by hand, so a numerical method is required

---

Step 4: Write a root-finding function

• Write a function that numerically solves: [ f(x) = 0]

• The function should:

  1. Start from an initial guess
  2. Iteratively update the estimate
  3. Stop when convergence is reached

---

Step 5: Output requirements

Your function must return:

• The estimated root of ( f(x) = 0 )
• The value of ( f(x) ) at the estimated root
• The number of iterations used

---

Workflow for Question (b)

Step 6: Series computation

• Identify the given series

• Write a function that:

  • Computes the sum of the series
  • Uses a specified number of terms or a stopping rule

• Return the final sum

---

One-line workflow summary

Define → Visualize → Differentiate → Solve numerically → Report results → Compute series sum

---

# 4.(a) Draw functions l(x) and f(x) on interval -4 < x < 1
l_func <- function(x) x^4 + 4*x^3 + 2*x^2 + 12*x + 5
f_func <- function(x) 4*x^3 + 12*x^2 + 4*x + 12
f_prime <- function(x) 12*x^2 + 24*x + 4

# Create sequence for plotting
x_vals <- seq(-4, 1, length.out = 1000)

# Plot both functions
plot(x_vals, l_func(x_vals), type = "l", col = "blue", 
     xlab = "x", ylab = "Function Value",
     main = "Functions l(x) and f(x)",
     ylim = range(c(l_func(x_vals), f_func(x_vals))))
lines(x_vals, f_func(x_vals), col = "red")

abline(h = 0, lty = 2, col = "gray")
legend("topright", legend = c("l(x)", "f(x)"), 
       col = c("blue", "red"), lty = 1)

cat("4.(a) Plot created.\n\n")

4.(a) Plot created.

# 4.(a)(ii) Function to solve transcendental equation using Newton's method
solve_transcendental <- function(x0 = 0, tol = 1e-10, max_iter = 1000) {
  x <- x0
  iterations <- 0
  
  for(i in 1:max_iter) {
    fx <- f_func(x)
    fpx <- f_prime(x)
    
    # Newton update: x_new = x - f(x)/f'(x)
    x_new <- x - fx / fpx
    
    # Check convergence
    if(abs(x_new - x) < tol && abs(fx) < tol) {
      break
    }
    
    x <- x_new
    iterations <- i
  }
  
  return(list(
    solution = x,
    function_value = f_func(x),
    iterations = iterations
  ))
}

# Find solution (try different starting points)
result <- solve_transcendental(x0 = -2)
cat("4.(a)(ii) Solution using Newton's method:\n")

4.(a)(ii) Solution using Newton's method:

cat("Estimated solution:", result$solution, "\n")

Estimated solution: -3 

cat("Function value at estimate:", result$function_value, "\n")

Function value at estimate: -6.583178e-12 

cat("Number of iterations:", result$iterations, "\n\n")

Number of iterations: 8 

# 4.(b) Function to sum series for log(1+x)
log_series <- function(x, t0 = 1e-10, max_terms = 1000) {
  # Check validity: -1 < x < 1
  if(x <= -1 || x >= 1) {
    stop("x must be in the interval (-1, 1)")
  }
  
  term <- x
  sum_series <- term
  n <- 1
  
  while(abs(term) >= t0 && n < max_terms) {
    n <- n + 1
    term <- (-1)^(n-1) * x^n / n
    sum_series <- sum_series + term
  }
  
  return(list(
    x = x,
    series_sum = sum_series,
    true_value = log(1 + x),
    terms_used = n,
    absolute_error = abs(sum_series - log(1 + x))
  ))
}

# Test the function
cat("4.(b) Testing log series function:\n")

4.(b) Testing log series function:

test_values <- c(0.2, -0.3, 0.9)

for(val in test_values) {
  result <- tryCatch({
    log_series(val)
  }, error = function(e) {
    list(x = val, error = e$message)
  })
  
  if("error" %in% names(result)) {
    cat("x =", val, ":", result$error, "\n")
  } else {
    cat("x =", val, ":\n")
    cat("Series sum:", result$series_sum, "\n")
    cat("True log(1+x):", result$true_value, "\n")
    cat("Terms used:", result$terms_used, "\n")
    cat("Absolute error:", result$absolute_error, "\n")
  }
}

x = 0.2 :
Series sum: 0.1823216 
True log(1+x): 0.1823216 
Terms used: 13 
Absolute error: 9.863138e-12 
x = -0.3 :
Series sum: -0.3566749 
True log(1+x): -0.3566749 
Terms used: 17 
Absolute error: 3.008227e-11 
x = 0.9 :
Series sum: 0.6418539 
True log(1+x): 0.6418539 
Terms used: 170 
Absolute error: 4.623013e-11 

5. Use iris dataset to compute:

1. means of the numeric columns
2. the means of the numeric columns by Species.

# (i) Means of numeric columns
numeric_cols <- iris[, 1:4]
col_means <- colMeans(numeric_cols)
cat("(i) Means of numeric columns:\n")

(i) Means of numeric columns:

cat("Sepal Length:", col_means[1], "\n")

Sepal Length: 5.843333 

cat("Sepal Width:", col_means[2], "\n")

Sepal Width: 3.057333 

cat("Petal Length:", col_means[3], "\n")

Petal Length: 3.758 

cat("Petal Width:", col_means[4], "\n\n")

Petal Width: 1.199333 

# (ii) Means by Species
cat("(ii) Means by Species:\n")

(ii) Means by Species:

species_means <- aggregate(. ~ Species, data = iris, FUN = mean)
print(species_means)

print(iris)

# 5. Iris dataset analysis
# (i) Overall means
cat("(i) Overall means of numeric columns:\n")

(i) Overall means of numeric columns:

cat("-------------------------------------\n")

-------------------------------------

overall <- colMeans(iris[1:4])
names(overall) <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
print(overall)

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 

# (ii) Means by Species
cat("(ii) Means by Species:\n")

(ii) Means by Species:

cat("----------------------\n")

----------------------

species <- c("setosa", "versicolor", "virginica")

for(sp in species) {
  cat(sp, ":\n")
  subset_means <- colMeans(iris[iris$Species == sp, 1:4]) #iris[rows, columns]
  print(subset_means)
}

setosa :
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       5.006        3.428        1.462        0.246 
versicolor :
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       5.936        2.770        4.260        1.326 
virginica :
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       6.588        2.974        5.552        2.026 

2023

1.(a) Make 3 vectors: (no loops)

Vector A: (3, 3.1, …, 6)

Vector B: (2, 2²/2, 2³/3, …, 2²⁵/25)

Vector C: Repeat (4,6,3) to get 10 total 4s Sum each vector.

#A
A <- seq(3,6,0.1)
A

 [1] 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6
[18] 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0

#B
n <- 1:25
B <- (2^n) / n
B

 [1] 2.000000e+00 2.000000e+00 2.666667e+00 4.000000e+00 6.400000e+00
 [6] 1.066667e+01 1.828571e+01 3.200000e+01 5.688889e+01 1.024000e+02
[11] 1.861818e+02 3.413333e+02 6.301538e+02 1.170286e+03 2.184533e+03
[16] 4.096000e+03 7.710118e+03 1.456356e+04 2.759411e+04 5.242880e+04
[21] 9.986438e+04 1.906502e+05 3.647221e+05 6.990507e+05 1.342177e+06

#C
C <- rep(c(4,6,3), 10)
C

 [1] 4 6 3 4 6 3 4 6 3 4 6 3 4 6 3 4 6 3 4 6 3 4 6 3 4 6 3 4 6 3

1.(b) Generate xVec, yVec with seed 50. Calculate: [No loops] xVec <- sample(0:999, 250, replace=T) yVec <- sample(0:999, 250, replace=T) x,y belongs to X,Y i) Sum(y − x) ii) Root over (Average of x²) iii) Sum of exp(−y_(i+1))/((x_i)+10) for i=1 to n−1

set.seed(50)
xVec <- sample(0:999, 250, replace=T)
yVec <- sample(0:999, 250, replace=T)
#(i)
ans_i <- sum(yVec-xVec)
ans_i

[1] -12273

#(ii)
ans_ii <- sqrt(mean(xVec^2))
ans_ii

[1] 600.7198

#(iii)
n <- length(xVec)
# Use 2:n for y (i+1) and 1:(n-1) for x (i)
ans_iii <- sum(exp(-yVec[2:n]) / (xVec[1:(n-1)] + 10))
ans_iii

[1] 0.004409429

2.(a) With seed 209: Generate X ~ Binomial(20, 0.5) and ε ~ N(0, 3), 100 each. Compute Y = 20 + 1.5X + ε. Save as x, y.

set.seed(209)
x <- rbinom(100, 20, 0.5)
ε <- rnorm(100, 0, 3)
y <- 20 + 1.5*x + ε

2.(b) Write function: input x, y → output matrix[a, b]ᵀ (col matrix…not transpose). b = (Σxy − (ΣxΣy)/n)/(Σx^2-((Σx)2)/n) a = ybar - b*xbar Store as mat1. & Compare to [20, 1.5]ᵀ which is the matrix of regression parameters of (a).

# --- Part 2.(b) ---
# Define the function
get_coefficients <- function(x, y) {
  n <- length(x)
  # Calculate necessary sums
  sum_x <- sum(x)
  sum_y <- sum(y)
  sum_xy <- sum(x * y)
  sum_x_sq <- sum(x^2)
  # Calculate slope b using the provided formula
  # Numerator: Σxy - (ΣxΣy)/n
  numerator <- sum_xy - (sum_x * sum_y) / n
  # Denominator: Σx^2 - (Σx)^2/n
  denominator <- sum_x_sq - (sum_x^2) / n
  b <- numerator / denominator
  # Calculate intercept a
  # a = ybar - b * xbar
  a <- mean(y) - b * mean(x)
  # Return as a column matrix (2 rows, 1 column)
  return(matrix(c(a, b), nrow = 2, ncol = 1))
}

# Store the result as mat1
mat1 <- get_coefficients(x, y)

# Define the theoretical parameters [20, 1.5] for comparison
theoretical <- matrix(c(20, 1.5), nrow = 2, ncol = 1)

# Display the results
print("Estimated Matrix (mat1):")

[1] "Estimated Matrix (mat1):"

print(mat1)

          [,1]
[1,] 20.695560
[2,]  1.457999

print("Theoretical Matrix:")

[1] "Theoretical Matrix:"

print(theoretical)

     [,1]
[1,] 20.0
[2,]  1.5

print("Difference (Estimated - Theoretical):")

[1] "Difference (Estimated - Theoretical):"

print(mat1 - theoretical)

            [,1]
[1,]  0.69556007
[2,] -0.04200109

2.(c) Create D = [1(column full of ones), x] and a column matrix M of 1000 rows using vector y. Compute (DᵀD)⁻¹DᵀM → store as mat2. Compare to [20, 1.5]ᵀ.

D <- cbind(1, x)  # First column = 1, second column = X
M <- matrix(rep(y, length.out = 1000), ncol = 1)

# Calculate (D^T D)^{-1} D^T M
# But D is 100x2, M is 1000x1 - dimensions don't match!
# This seems like an error. Let me assume M should be based on y (100x1)
M_correct <- matrix(y, ncol = 1)  # 100x1 matrix

mat2 <- solve(t(D) %*% D) %*% t(D) %*% M_correct
cat("Estimated parameters (mat2):\n")

Estimated parameters (mat2):

print(mat2)

       [,1]
  20.695560
x  1.457999

cat("Difference from true parameters:\n")

Difference from true parameters:

print(mat2 - matrix(c(20, 1.5), nrow = 2, ncol = 1))

         [,1]
   0.69556007
x -0.04200109

2.(d) X <- matrix(c(0.63, 0.75, 1.53, -1.70, -0.93, 1.23, 0.80, 1.11, 2.84), nrow = 3, byrow = TRUE) Replace: Row 2 with (0.47, 2.33, 87) Row1Col3=63, Row3Col3=75 Find column means and SDs.

x <- matrix(c(0.63, 0.75, 1.53,
                    -1.70, -0.93, 1.23,
                     0.80, 1.11, 2.84), nrow = 3, byrow = TRUE)
x[2,] <- c(0.47, 2.33, 87)
x[1,3] = 63
x[3,3] = 75
x

     [,1] [,2] [,3]
[1,] 0.63 0.75   63
[2,] 0.47 2.33   87
[3,] 0.80 1.11   75

col_means <- apply(x, 2, mean)
col_means

[1]  0.6333333  1.3966667 75.0000000

col_sds <- apply(x, 2, sd)
col_sds

[1]  0.1650253  0.8280902 12.0000000

1. X ~ Poisson(λ=2.7). Find:

1. P(X > 3)
2. P(2 < X < 8)
3. Smallest k: P(X ≤ k) ≥ 0.5

#(i)
px3_ <- 1-ppois(3,2.7)
px3_

[1] 0.2859078

#(ii)
px2_8 <- sum(dpois(3:7,2.7))
px2_8

[1] 0.4997543

#(iii)
k <- qpois(0.5, 2.7)
k

[1] 3

2. Generate 3000 obs from Poisson(λ=1.5) with seed 101. Estimate:

1. P(X > 3)
2. P(2 < X < 8)
3. p(x<=3)

set.seed(101)
sim_x <- rpois(3000, 1.5)

# (i) Estimate P(X > 3)
prob_i <- mean(sim_x > 3)
print(paste("Estimate for P(X > 3):", prob_i))

[1] "Estimate for P(X > 3): 0.071"

# (ii) Estimate P(2 < X < 8)
prob_ii <- mean(sim_x > 2 & sim_x < 8)
print(paste("Estimate for P(2 < X < 8):", prob_ii))

[1] "Estimate for P(2 < X < 8): 0.193666666666667"

# (iii) Estimate P(X <= 3)
prob_iii <- mean(sim_x <= 3)
print(paste("Estimate for P(X <= 3):", prob_iii))

[1] "Estimate for P(X <= 3): 0.929"

3. X ~ Exponential(mean=1). Find:

1. P(0.5 < X < 1.5)
2. k such that P(X < k) = 0.8

px0.5_1.5 <- pexp(1.5, rate=1)-pexp(0.5, rate=1) #rate = 1/mean = 1/lambda
px0.5_1.5

[1] 0.3834005

k <- qexp(0.8, rate=1)
k

[1] 1.609438

4. Generate 5000 obs from Exp(1) with seed 102. Estimate:

1. P(0.5 < X < 1.5)
2. P(X > 1.61)

set.seed(102)
# Note: Since mean = 1, rate = 1/mean = 1
sim_x_exp <- rexp(5000, rate = 1)
# (i) Estimate P(0.5 < X < 1.5)
prob_i_exp <- mean(sim_x_exp > 0.5 & sim_x_exp < 1.5)
print(paste("Estimate for P(0.5 < X < 1.5):", prob_i_exp))

[1] "Estimate for P(0.5 < X < 1.5): 0.3832"

# (ii) Estimate P(X > 1.61)
prob_ii_exp <- mean(sim_x_exp > 1.61)
print(paste("Estimate for P(X > 1.61):", prob_ii_exp))

[1] "Estimate for P(X > 1.61): 0.1924"

1. Write function: Add two vectors by truncating longer to match shorter length.

trunc_add <- function(v1, v2) {
  n <- min(length(v1), length(v2))
  return(v1[1:n] + v2[1:n])
}

2. Write A.sort (sorts ascending) and own.med (finds median using A.sort). With seed 200, generate 500 obs from N(10,10), find median using own.med.

# 1. Define Functions
A.sort <- function(x) {
  return(sort(x)) # Uses R's built-in sort (ascending)
}

own.med <- function(x) {
  sorted_x <- A.sort(x)
  n <- length(sorted_x)
  if (n %% 2 == 1) {
    return(sorted_x[(n + 1) / 2])
  } else {
    return(mean(sorted_x[c(n / 2, n / 2 + 1)]))
  }
}

# 2. Generate Data and Test
set.seed(200)
# Assuming N(10, 10) means mean=10, sd=10 based on previous context
data <- rnorm(500, mean = 10, sd = 10) 
result <- own.med(data)
print(paste("Median:", result))

[1] "Median: 10.1536202073796"

3. Plot f(x)=x³−2x²−4. Use Newton-Raphson to solve f(x)=0. Function inputs: x0, tolerance t0. Outputs: root, f(root), iteration count.

# 1. Define the function f(x) and its derivative df(x)
f <- function(x) x^3 - 2*x^2 - 4
df <- function(x) 3*x^2 - 4*x

# 2. Plot f(x)
curve(f, from = -2, to = 4, main = "Plot of f(x)")
abline(h = 0, col = "red") # Shows the x-axis

# 3. Newton-Raphson Function
newton_solve <- function(x0, t0) {
  iter <- 0
  x <- x0
  while (abs(f(x)) > t0) {
    x <- x - f(x) / df(x)
    iter <- iter + 1
  }
  return(list(root = x, value = f(x), iterations = iter))
}

# Example Usage (Guessing 2)
newton_solve(2, 0.0001)

$root
[1] 2.594318

$value
[1] 4.882565e-05

$iterations
[1] 4

1. In birthwt dataset: treat low, race, smoke, ht, ui as categorical. Find mean of numeric columns.

# Load data
library(MASS)
data(birthwt)

birthwt$low <- factor(birthwt$low) #Factor converts numeric/text data into categories (groups).
birthwt$race <- factor(birthwt$race)
birthwt$smoke <- factor(birthwt$smoke)
birthwt$ht <- factor(birthwt$ht)
birthwt$ui <- factor(birthwt$ui)

numeric_cols <- c("age", "lwt", "bwt", "ptl", "ftv")
numeric_means <- colMeans(birthwt[numeric_cols])

cat("(a) Means of numeric columns:\n")

(a) Means of numeric columns:

print(numeric_means)

         age          lwt          bwt          ptl          ftv 
  23.2380952  129.8148148 2944.5873016    0.1957672    0.7936508 

2. Write function: For bwt variable, compute mean and SD for each race group.

bwt_stats_by_race <- function(data) {
  races <- levels(data$race)
  results <- data.frame()
  
  for(r in races) {
    bwt_values <- data$bwt[data$race == r]
    mean_val <- mean(bwt_values)
    sd_val <- sd(bwt_values)
    results <- rbind(results, data.frame(Race = r, Mean = mean_val, SD = sd_val))
  }
  
  return(results)
}

cat("\n(b) BWT statistics by race:\n")

(b) BWT statistics by race:

print(bwt_stats_by_race(birthwt))