Tugas Pra UTS APG Kelompok 5

Membuat fungsi dari metode-metode yang telah dipelajari (pertemuan 1-6). Fungsi diuji kelayakannya dengan menggunakan dataset yang tersedia dalam r.

Anggota

• Alif Andika Putra, 16.8994
• M Ghozy Al Haqqoni, 16.9293
• M Rifqi Mubarak, 16.9304
• Mukhammad Kharis, 16.9311
• Novia Permatasari, 16.9335

R Dataset

Dataset yang digunakan untuk menguji kelayakan yaitu dataset state.x77. Berikut contoh tampilan dataset:

dataset <- state.x77
head(dataset)

           Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
California      21198   5114        1.1    71.71   10.3    62.6    20 156361
Colorado         2541   4884        0.7    72.06    6.8    63.9   166 103766

Mean Vector

#Fungsi
vectorMean <- function(x) {
  n = dim(x)[1]
  xbar = rep(1,n) %*% x / n
  return(xbar)
}

Berikut penggunaan fungsi yang telah dibuat dengan dataset.

vectorMean(dataset)

     Population Income Illiteracy Life Exp Murder HS Grad  Frost     Area
[1,]    4246.42 4435.8       1.17  70.8786  7.378  53.108 104.46 70735.88

Selanjutnya dilakukan pengujian terhadap fungsi colMeans yang telah tersedia dalam R. Fungsi valid jika vectorMean(dataset) == colMeans(dataset).

#Validitas
if(colMeans(dataset) == vectorMean(dataset)){
  print("Fungsi valid")
}else
  print("Fungsi tidak valid")

[1] "Fungsi valid"

Covariance and Correlation Matrix

Covariance Matrix

#Fungsi
matrixCov <- function(x) {
  n = dim(x)[1]
  ones = rep(1,n)
  c = diag(1,n) - ones %*% t(ones) / n
  xc = c %*% x
  matCov = t(xc) %*% xc / (n-1)
  return(matCov)
}

Berikut penggunaan fungsi yang telah dibuat dengan dataset.

matrixCov(dataset)

              Population        Income   Illiteracy      Life Exp       Murder       HS Grad
Population 19931683.7588   571229.7796  292.8679592 -4.078425e+02  5663.523714  -3551.509551
Income       571229.7796   377573.3061 -163.7020408  2.806632e+02  -521.894286   3076.768980
Illiteracy      292.8680     -163.7020    0.3715306 -4.815122e-01     1.581776     -3.235469
Life Exp       -407.8425      280.6632   -0.4815122  1.802020e+00    -3.869480      6.312685
Murder         5663.5237     -521.8943    1.5817755 -3.869480e+00    13.627465    -14.549616
HS Grad       -3551.5096     3076.7690   -3.2354694  6.312685e+00   -14.549616     65.237894
Frost        -77081.9727     7227.6041  -21.2900000  1.828678e+01  -103.406000    153.992163
Area        8587916.9494 19049013.7510 4018.3371429 -1.229410e+04 71940.429959 229873.192816
                  Frost          Area
Population -77081.97265  8.587917e+06
Income       7227.60408  1.904901e+07
Illiteracy    -21.29000  4.018337e+03
Life Exp       18.28678 -1.229410e+04
Murder       -103.40600  7.194043e+04
HS Grad       153.99216  2.298732e+05
Frost        2702.00857  2.627039e+05
Area       262703.89306  7.280748e+09

Selanjutnya dilakukan pengujian terhadap fungsi cov yang telah tersedia dalam R. Fungsi valid jika matrixCov(dataset) == cov(dataset).

#Validitas
if(cov(dataset) == matrixCov(dataset)){
  print("Fungsi valid")
}else
  print("Fungsi tidak valid")

[1] "Fungsi valid"

Correlation Matrix

#Fungsi
matrixCor <- function(x) {
  matCov = matrixCov(x)
  result = matrix(nrow = dim(matCov)[1], ncol = dim(matCov)[2])
  for (i in 1:dim(matCov)[1]) {
    for (j in 1:dim(matCov)[2]) {
      result[i,j] = matCov[i,j] / (sqrt(matCov[i,i]) * sqrt(matCov[j,j]))
    }
  }
  return(result)
}

Berikut menguji kelayakan dengan dataset. Set Variable yang digunakan untuk melihat matriks korelasi yaitu variabel Population, Income, serta Illiteracy.

matrixCor(dataset[,1:3])

          [,1]       [,2]       [,3]
[1,] 1.0000000  0.2082276  0.1076224
[2,] 0.2082276  1.0000000 -0.4370752
[3,] 0.1076224 -0.4370752  1.0000000

#Cek validitas fungsi dengan fungsi cor()
cor(dataset[,1:3])

           Population     Income Illiteracy
Population  1.0000000  0.2082276  0.1076224
Income      0.2082276  1.0000000 -0.4370752
Illiteracy  0.1076224 -0.4370752  1.0000000

Contour Plot

Fungsi yang dibuat menunjukkan garis kontur (plot).

require(ellipse)
#Fungsi
kontur <- function(data){
  plot(data, main = "Contour Plot")
  rata<-vectorMean(data)
  cov<-matrixCov(data)
  lines( ellipse( cov, centre = rata,level=0.20) , col='black')
  lines( ellipse( cov, centre = rata,level=0.50) , col='black')
  lines( ellipse( cov, centre = rata,level=0.75) , col='blue')
  lines( ellipse( cov, centre = rata,level=0.90) , col='blue')
  lines( ellipse( cov, centre = rata,level=0.95) , col='red')
  lines( ellipse( cov, centre = rata,level=0.99) , col='red')
}

Berikut menguji kelayakan dengan dataset.

par(mfrow = c(2,2))
for(i in 1:4){
kontur(dataset[,i:5])
}
par(mfrow = c(1,1))

T2-Hotelling

Satu Sampel

#Fungsi
tHotellingTest1Sampel <- function(x, miu, alpha) {
  n = dim(x)[1]
  p = dim(x)[2]
  ybar = t(vectorMean(x))
  s = matrixCov(x)
  t2 = n * t(ybar - miu) %*% solve(s) %*% (ybar - miu)
  t2table = (n - 1) * p / ((n - 1) - p + 1) * qf(1 - alpha, p, (n - 1) - p + 1)
  rejectH0 = F
  if (t2 >= t2table) {
    rejectH0 = T
  }
  #output
  writeLines('')
  writeLines('\tUji T2 Hotelling Satu Sampel')
  writeLines('')
  writeLines(paste('T2 Obs =', t2))
  writeLines(paste('T2 Tab =', t2table))
  if (rejectH0) {
    writeLines('Keputusan = Tolak H0')
  } else {
    writeLines('Keputusan = Gagal Tolak H0')
  }
}

Berikut menguji kelayakan dengan dataset.

miu <- colMeans(dataset) #Karena colMeans == miu -> keputusan Gagal Tolak H0
tHotellingTest1Sampel(dataset, miu, .05)

    Uji T2 Hotelling Satu Sampel

T2 Obs = 1.54269109497771e-27
T2 Tab = 20.2357556026252
Keputusan = Gagal Tolak H0

Dua Sampel Dependen

#Fungsi
tHotellingTest2SampelDep <- function(x, y, alpha) {
  d = x - y
  n = dim(d)[1]
  p = dim(d)[2]
  dbar = t(vectorMean(d))
  s = matrixCov(d)
  t2 = n * t(dbar) %*% solve(s) %*% (dbar)
  t2table = (n - 1) * p / ((n - 1) - p + 1) * qf(1 - alpha, p, (n - 1) - p + 1)
  rejectH0 = F
  if (t2 >= t2table) {
    rejectH0 = T
  }
  #output
  writeLines('')
  writeLines('\tUji T2 Hotelling Dua Sampel Dependen')
  writeLines('')
  writeLines(paste('T2 Obs =', t2))
  writeLines(paste('T2 Tab =', t2table))
  if (rejectH0) {
    writeLines('Keputusan = Tolak H0')
  } else {
    writeLines('Keputusan = Gagal Tolak H0')
  }
}

Berikut menguji kelayakan dengan dataset.

miu <- colMeans(dataset)
a1 <- dataset[1:25,1:3]
a2 <- dataset[26:50,1:3]
tHotellingTest2SampelDep(a1, a2, .05)

    Uji T2 Hotelling Dua Sampel Dependen

T2 Obs = 4.02788412929234
T2 Tab = 9.97895450831698
Keputusan = Gagal Tolak H0

Dua Sampel Independen

#Fungsi
tHotellingTest2SampelInd <- function(x, y, alpha) {
  nx = dim(x)[1]
  px = dim(x)[2]
  ny = dim(y)[1]
  py = dim(y)[2]
  xbar = t(vectorMean(x))
  ybar = t(vectorMean(y))
  sx = matrixCov(x)
  sy = matrixCov(y)
  sp = ((nx - 1) * sx + (ny - 1) * sy) / (nx + ny - 2)
  t2 = nx * ny / (nx + ny) * t(xbar - ybar) %*% solve(sp) %*% (xbar - ybar)
  t2table = (nx + ny - 2) * px / ((nx + ny - 2) - px + 1) * qf(1 - alpha, px, (nx + ny - 2) - px + 1)
  rejectH0 = F
  if (t2 >= t2table) {
    rejectH0 = T
  }
  #output
  writeLines('')
  writeLines('\tUji T2 Hotelling Dua Sampel Independen')
  writeLines('')
  writeLines(paste('T2 Obs =', t2))
  writeLines(paste('T2 Tab =', t2table))
  if (rejectH0) {
    writeLines('Keputusan = Tolak H0')
  } else {
    writeLines('Keputusan = Gagal Tolak H0')
  }
}

Berikut menguji kelayakan dengan dataset.

tHotellingTest2SampelInd(a1, a2, .05)

    Uji T2 Hotelling Dua Sampel Independen

T2 Obs = 3.0625002928296
T2 Tab = 8.78664499452393
Keputusan = Gagal Tolak H0

Likelihood Ratio Test

Fungsi uji rasio likelihood digunakan untuk dataset satu sampel.

#Fungsi
likelihoodTest <- function(x , xbar_0, alfa){
  p <- ncol(x)
  n <- nrow(x)
  v <- n - 1
  
  xbar <- vectorMean(x)
  
  sgm <- matrixCov(x)
  sgm0 <- matrix(rep(0, (p^2)) , p)    
  for (i in 1:p) {
    for (k in 1:p) {
      sum <- 0
      for (j in 1:n) {
        sum <- sum + ((x[j, i] - xbar_0[i]) * (x[j, k] - xbar_0[k]))
      }
      sgm0[i, k] <- sum / (n - 1)
    }
  }
  sgm0
  
  # wilk <- mpower(det(sgm)/det(sgm0), n/2)
  t2table <- qf((1-alfa) , p, v-p+1) * (v*p) / (v-p+1)
  t2 <- (n-1)*(det(sgm0)/det(sgm)) - (n-1)
  
  rejectH0 = F
  if(t2 >= t2table){
    rejectH0 = T
  }
  #output
  writeLines('')
  writeLines('\tUji Rasio Likelihood')
  writeLines('')
  writeLines(paste('T2 Obs =', t2))
  writeLines(paste('T2 Tab =', t2table))
  if (rejectH0) {
    writeLines('Keputusan = Tolak H0')
  } else {
    writeLines('Keputusan = Gagal Tolak H0')
  }
}

Berikut menguji kelayakan dengan dataset. Fungsi valid jika Keputusan sama dengan keputusan menggunakan fungsi T2-Hotelling Satu Sampel yaitu Gagal Tolak H0.

likelihoodTest(dataset, miu, .05)

    Uji Rasio Likelihood

T2 Obs = 0
T2 Tab = 20.2357556026252
Keputusan = Gagal Tolak H0

TERIMAKASIH. Mohon maaf atas segala kekurangan.