Nomor 1
# Input Parameter
mu <- c(0, 2)
Sigma <- matrix(c(2, sqrt(2)/2, sqrt(2)/2, 1), 2, 2)
# Generate Random Bivariate Normal
library(MASS)
bvn1 <- mvrnorm(500, mu = mu, Sigma = Sigma )
head(bvn1, 10)
## [,1] [,2]
## [1,] -1.03950381 1.441832
## [2,] -0.90705563 1.164982
## [3,] 3.18226021 2.802930
## [4,] -1.02797812 2.592629
## [5,] 2.79600094 2.706602
## [6,] -1.49752027 3.085969
## [7,] 0.84078650 1.924910
## [8,] 0.58283599 1.694129
## [9,] 0.35528058 2.648239
## [10,] -0.06721738 1.873636
# (b) Ekspresi jarak kuadrat
jarak_kuadrat <- mahalanobis(bvn1, center = colMeans(bvn1), cov = cov(bvn1))
head(jarak_kuadrat,10)
## [1] 0.64758708 0.82321127 5.44476365 1.72017514 4.19820047 4.52424753
## [7] 0.55143728 0.52786857 0.40138462 0.02161567
# (c) Kontur ellips dengan probabilitas 50%
library(mixtools)
plot(bvn1, xlab="X1",ylab="X2", main = "Kontur Ellipsoid Normal Bivariat
dengan Probabilitas 50%")
ellipse(mu, Sigma, 0.5, col="red", lwd=2)
Nomor 2
# import data
library(readxl)
PDRB_ADHB_2023 <- read_excel("C:/Users/anggi/Documents/KULIAH/SEMESTER 5/APG/PRAKTIKUM 3/PDRB_ADHB_2023.xlsx")
x1 <- PDRB_ADHB_2023$P
x2 <- PDRB_ADHB_2023$Q
x3 <- PDRB_ADHB_2023$R
pdrb <- cbind(x1,x2,x3)
# Cek Normalitas Data Multivariat
library(MVN)
#uji Royston
UjiRoy <- mvn(pdrb, mvnTest = "royston", multivariatePlot = "qq")
UjiRoy
## $multivariateNormality
## Test H p value MVN
## 1 Royston 57.03905 1.042193e-13 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling x1 6.3353 <0.001 NO
## 2 Anderson-Darling x2 5.8808 <0.001 NO
## 3 Anderson-Darling x3 8.7529 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
## x1 38 17388.079 29799.90 7988.0 656 150213 3025.25 13921.75 2.917140 8.746808
## x2 38 7121.500 12713.33 3453.5 396 76588 1439.50 6957.00 4.327410 20.553914
## x3 38 9915.474 26027.01 2528.0 133 149694 956.75 5206.50 4.274544 19.275075
#Uji Mardia
UjiMar <- mvn(pdrb, mvnTest = "mardia", multivariatePlot = "qq")
UjiMar
## $multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 263.897528346024 6.45735543090546e-51 NO
## 2 Mardia Kurtosis 22.4785207610662 0 NO
## 3 MVN <NA> <NA> NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling x1 6.3353 <0.001 NO
## 2 Anderson-Darling x2 5.8808 <0.001 NO
## 3 Anderson-Darling x3 8.7529 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
## x1 38 17388.079 29799.90 7988.0 656 150213 3025.25 13921.75 2.917140 8.746808
## x2 38 7121.500 12713.33 3453.5 396 76588 1439.50 6957.00 4.327410 20.553914
## x3 38 9915.474 26027.01 2528.0 133 149694 956.75 5206.50 4.274544 19.275075
# Optional: Transformasi jika data tidak normal
library(car)
pow <- powerTransform (pdrb)
summary(pow)
## bcPower Transformations to Multinormality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## x1 0.0778 0 -0.0660 0.2217
## x2 0.0391 0 -0.1160 0.1942
## x3 -0.0080 0 -0.1187 0.1027
##
## Likelihood ratio test that transformation parameters are equal to 0
## (all log transformations)
## LRT df pval
## LR test, lambda = (0 0 0) 2.800423 3 0.42343
##
## Likelihood ratio test that no transformations are needed
## LRT df pval
## LR test, lambda = (1 1 1) 299.6812 3 < 2.22e-16
pdrb_t <- log(pdrb)
#Uji Royston
UjiRoy_t <- mvn(pdrb_t, mvnTest = "royston", multivariatePlot = "qq")
UjiRoy_t
## $multivariateNormality
## Test H p value MVN
## 1 Royston 0.5131584 0.6204803 YES
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling x1 0.3056 0.5505 YES
## 2 Anderson-Darling x2 0.1604 0.9435 YES
## 3 Anderson-Darling x3 0.3562 0.4397 YES
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew
## x1 38 8.877243 1.332537 8.985427 6.486161 11.91981 8.014525 9.541152 0.1879969
## x2 38 8.156595 1.162994 8.146830 5.981414 11.24620 7.271660 8.847470 0.2354952
## x3 38 7.842287 1.532502 7.829989 4.890349 11.91635 6.863488 8.557288 0.4682832
## Kurtosis
## x1 -0.4082591
## x2 -0.2332129
## x3 0.1489837
Nomor 3
# Data perusahaan terbesar
data <- matrix(c(108.28, 17.05, 1484.10,
152.36, 16.59, 750.33,
95.04, 10.91, 766.42,
65.45, 14.14, 1110.46,
62.97, 9.52, 1031.29,
263.99, 25.33, 195.26,
265.19, 18.54, 193.83,
285.06, 15.73, 191.11,
92.01, 8.10, 1175.16,
165.68, 11.13, 211.15), ncol = 3, byrow = TRUE)
colnames(data) <- c("sales", "profits", "assets")
# (a) Pemeriksaan normalitas multivariat
# Uji Roytoson
library(MVN)
result <- mvn(data, mvnTest = "royston", multivariatePlot = "qq")
result
## $multivariateNormality
## Test H p value MVN
## 1 Royston 3.197212 0.1437248 YES
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling sales 0.5944 0.0887 YES
## 2 Anderson-Darling profits 0.2713 0.5890 YES
## 3 Anderson-Darling assets 0.5669 0.1053 YES
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## sales 10 155.603 86.466486 130.320 62.97 285.06 92.7675 239.4125
## profits 10 14.704 5.117647 14.935 8.10 25.33 10.9650 16.9350
## assets 10 710.911 486.882193 758.375 191.11 1484.10 199.2325 1090.6675
## Skew Kurtosis
## sales 0.4094674 -1.6803972
## profits 0.5589459 -0.6667657
## assets 0.1027204 -1.7232822
# Uji Shapiro
library(RVAideMemoire)
mshapiro.test(data)
##
## Multivariate Shapiro-Wilk normality test
##
## data: (sales,profits,assets)
## W = 0.92772, p-value = 0.4258
# (b) Menghitung jarak Mahalanobis untuk perusahaan baru
new_company <- c(230, 24, 230.5)
mean_vector <- colMeans(data)
cov_matrix <- cov(data)
# Jarak Mahalanobis
mahalanobis_distance <- mahalanobis(new_company, center = mean_vector, cov = cov_matrix)
# Periksa apakah berada di dalam atau di luar kontur 95%
if(mahalanobis_distance <= qchisq(0.95, df = 3)) {
print("Perusahaan berada di dalam kontur 95%.")
} else {
print("Perusahaan berada di luar kontur 95%.")
}
## [1] "Perusahaan berada di dalam kontur 95%."