regresi_berganda <- function(X, y) {
  if (nrow(X) != length(y)) {
    stop("Jumlah baris X harus sama dengan panjang y")
  }

  # X'X
  XtX <- t(X) %*% X
  
  # Invers dari X'X
  XtX_inv <- solve(XtX)
  
  # X'y
  Xty <- t(X) %*% y
  
  # Beta_hat: (X'X)^-1 X'y
  beta_hat <- XtX_inv %*% Xty
  return(beta_hat)
}

# Dataset
X <- matrix(c(1, 1, 0,
              1, 2, 1,
              1, 3, 1,
              1, 4, 0), ncol=3)
y <- c(5, 7, 8, 9)

# Regresi berganda
regresi_berganda(X, y)
##          [,1]
## [1,] 2.284553
## [2,] 1.869919
## [3,] 1.463415