資料集探索與分析

1
2
3
4
5

data = read.table('bank2.dat')
colnames(data) = c('X1','X2','X3','X4','X5','X6')
data = data.frame(data, genuine = c(rep(1, 100), rep(0, 100)))
attach(data)
head(data)

1
2
3
4
5

print('Quartiles:')
summary(data)

print('NAs:')
sapply(data, function(x) sum(is.na(x)))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

library(ggcorrplot)

cor_matrix = cor(data)

ggcorrplot(cor_matrix, 
           method = "circle",
           type = "full",
           lab = TRUE,
           lab_size = 3,
           colors = c("blue", "white", "red"),
           outline.color = "gray",
           legend.title = "Correlation",
           show.legend = TRUE)

1

pairs(data, upper.panel = NULL)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

pairs_plot = function(x, y, x_name, y_name){
  par(mar = c(5, 5, 4, 10))
  plot(x, y, main = "Swiss bank notes", xlab = x_name, ylab = y_name)
  points(x[1:100], y[1:100], col = "blue")
  points(x[101:200], y[101:200], col = "red")
  legend("topright", legend = c("Genuine", "Counterfeit"), col = c("blue", "red"), pch = 1, xpd = TRUE, inset = c(-0.3, 0))
}

num = combn(1:6, 2)
for (i in 1:ncol(num)) {
  x = paste0("X", num[1, ][i])
  y = paste0("X", num[2, ][i])
  pairs_plot(get(x), get(y), x, y)
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

library(lattice)
library(grid)

groups = data[, 7]

pch_types = c(1, 1)              
colors = c("blue", "red")        
pch_group = pch_types[ifelse(groups == 0, 1, 2)]
col_group = colors[ifelse(groups == 0, 1, 2)]

combos = combn(1:6, 3)

n_rows = 1        
n_cols = 1        
plots_per_page = n_rows * n_cols   

for (page in seq(1, ncol(combos), by = plots_per_page)) {
  grid.newpage()

  for (i in 0:(plots_per_page - 1)) {
    index = page + i
    if (index > ncol(combos)) break

    x_name = paste0("X", combos[1, index])
    y_name = paste0("X", combos[2, index])
    z_name = paste0("X", combos[3, index])

    x_data = data[[x_name]]
    y_data = data[[y_name]]
    z_data = data[[z_name]]

    p = cloud(z_data ~ x_data * y_data,
          pch = pch_group,
          col = col_group,
          cex = 1.2,
          ticktype = "detailed",
          main = paste("Swiss bank notes:", x_name, y_name, z_name),
          screen = list(z = -90, x = -90, y = 45),
          scales = list(arrows = FALSE, col = "black", distance = 1, cex = 0.5),
          xlab = list(x_name, rot = -10, cex = 1.2),
          ylab = list(y_name, rot = 10, cex = 1.2),
          zlab = list(z_name, rot = 90, cex = 1.1),

          par.settings = list(
            axis.line = list(col = "black"),    
            box.3d = list(
              col = c(
                "black", "transparent", "transparent","black", "black", "transparent", "transparent", "transparent",   
                "transparent", "transparent", "transparent", "transparent"))),

          key = list(
            space = "right",
            points = list(pch = pch_types, col = colors, cex = 1.5),
            text = list(c("Genuine", "Counterfeit")),
            border = FALSE
          )
        )

    print(p, split = c((i %% n_cols) + 1, (i %/% n_cols) + 1, n_cols, n_rows), more = TRUE)
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

library(plotly)
library(htmlwidgets)

x_data = data$X1
y_data = data$X2
z_data = data$X6
groups = data[, 7]

groups = factor(groups, levels = c(0, 1), labels = c("Genuine", "Counterfeit"))

x_range <- range(x_data)
y_range <- range(y_data)
z_range <- range(z_data)

x_scaled = (x_data - x_range[1]) / (x_range[2] - x_range[1])
y_scaled = (y_data - y_range[1]) / (y_range[2] - y_range[1])
z_scaled = (z_data - z_range[1]) / (z_range[2] - z_range[1])

tick_positions = seq(0, 1, length.out = 5)
tick_labels_x = round(seq(x_range[1], x_range[2], length.out = 5), 1)
tick_labels_y = round(seq(y_range[1], y_range[2], length.out = 5), 1)
tick_labels_z = round(seq(z_range[1], z_range[2], length.out = 5), 1)

hover_text <- paste0(
  "X1: ", round(x_data, 1), "<br>",
  "X2: ", round(y_data, 1), "<br>",
  "X6: ", round(z_data, 1)
)

fig <- plot_ly(
  x = x_scaled,
  y = y_scaled,
  z = z_scaled,
  color = groups,
  colors = c("blue", "red"),
  type = "scatter3d",
  mode = "markers",
  marker = list(size = 5),
  text = hover_text,
  hoverinfo = "text" 
)

fig <- fig %>% layout(
  title = "3D Scatter Plot of X1, X2, X6",
  scene = list(
    xaxis = list(title = "X1", range = c(0, 1), tickvals = tick_positions, ticktext = tick_labels_x),
    yaxis = list(title = "X2", range = c(0, 1), tickvals = tick_positions, ticktext = tick_labels_y),
    zaxis = list(title = "X6", range = c(0, 1), tickvals = tick_positions, ticktext = tick_labels_z)
  )
)

saveWidget(fig, "3d_scatter_plot.html")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

name = NULL
for (i in 1:6) {
  name = c(name, paste0("X", i, "_Genuine"))
  name = c(name, paste0("X", i, "_Counterfeit"))
}

groups = list()
for (i in 1:length(name)){
  if (i %% 2 == 1) {
    groups[[name[i]]] = data[1:100, as.integer(i / 2) + 1]
  } else {
    groups[[name[i]]] = data[101:200, as.integer(i / 2)]
  }
}

means <- sapply(groups, mean)

par(mfrow=c(1,3))
for (i in seq(1, length(name) ,by = 2)){
  boxplot(groups[i:(i + 1)], names = name[i:(i + 1)], frame = TRUE, main = "Swiss bank notes", cex.axis=0.99)
  
  for (j in 0:1) {
    lines(c(j + 0.6, j + 1.4), rep(means[[name[(i + j)]]], 2), lty = "dotted", lwd = 1.2)
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

data_Genuine = data[1:100, 1:6]   
data_Counterfeit = data[101:200, 1:6]

for (i in 1:6) {
  breaks_seq = seq(min(data[, i]), max(data[, i]), by = 0.1)
  min_x = floor(min(data[, i]))
  max_x = ceiling(max(data[, i]))

  hist(data_Genuine[, i],
       breaks = breaks_seq,
       col = rgb(0, 0, 1, 0.5),   
       border = "blue",
       xlim = c(min_x, max_x),
       ylim = c(0, max(table(cut(data[, i], breaks_seq))) + 5),
       main = "Histogram",
       xlab = paste0("X", i),
       axes = FALSE) 
  
  hist(data_Counterfeit[, i],
       breaks = breaks_seq,
       col = rgb(1, 0, 0, 0.5),   
       border = "red",
       add = TRUE)
  
  axis(side = 1, at = seq(min_x, max_x, by = 1))
  axis(side = 2)
  
  legend("topright",
         legend = c("Genuine", "Counterfeit", "Overlap"),
         fill = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5), rgb(0.6, 0, 0.4, 0.7)),
         border = c("blue", "red", "purple"))
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

library(KernSmooth)

for (i in 1:6) {
  
  data_Genuine = data[1:100, i]   
  data_Counterfeit = data[101:200, i]  

  fh1 = bkde(data_Genuine, kernel = "biweight")  
  fh2 = bkde(data_Counterfeit, kernel = "biweight")  

  x_min = floor(min(c(fh1$x, fh2$x)))
  x_max = ceiling(max(c(fh1$x, fh2$x)))
  y_max = max(c(fh1$y, fh2$y)) * 1.1

  x_ticks <- seq(x_min, x_max, by = 1)

  plot(fh1, type = "l", lwd = 2,
       xlab = "Counterfeit / Genuine", 
       ylab = paste0("Density estimates for X", i, sep = ""),
       col = "blue", main = paste("Swiss bank notes - X", i),
       xlim = c(x_min, x_max), ylim = c(0, y_max),
       axes = FALSE) 
  
  lines(fh2, lty = "dotted", lwd = 2, col = "red")

  axis(side = 1, at = x_ticks, labels = x_ticks)
  axis(side = 2)
  box()

  legend("topright",
         legend = c("Genuine", "Counterfeit"),
         col = c("blue", "red"),
         lty = c("solid", "dotted"),
         lwd = 2,
         cex = 0.8)
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

library(tourr)

x = data[1:200, ]
y = NULL

for (i in 1:6) {
  z = (x[, i] - min(x[, i])) / (max(x[, i]) - min(x[, i]))  # zero-one scaling
  y = cbind(y, z)
}

Type = data[, 7]
f = as.integer(Type)
grid = seq(0, 2 * pi, length = 1000)

plot(grid, 
     andrews(y[1, ])(grid), 
     type = "l", 
     lwd = 1.2, 
     main = "Andrews' curves (All Bank data)", 
     axes = FALSE, 
     frame = TRUE, 
     ylim = c(-0.5, 0.6), 
     ylab = "", 
     xlab = "", 
     col = ifelse(f[1] == 1, "blue", "red"), 
     lty = ifelse(f[1] == 1, "solid", "dotted"))

for (i in 2:200) {
  lines(grid, 
        andrews(y[i, ])(grid), 
        col = ifelse(f[i] == 1, "blue", "red"), 
        lwd = 1.2, 
        lty = ifelse(f[i] == 1, "solid", "dotted"))
}

axis(side = 2, at = seq(-0.5, 0.6, 0.2), labels = seq(-0.5, 0.6, 0.2))
axis(side = 1, at = seq(0, 7, 1), labels = seq(0, 7, 1))

legend("topright", 
       legend = c("Genuine", "Counterfeit"), 
       col = c("blue", "red"), 
       lwd = 1.5, 
       lty = c("solid", "dotted"))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# summary
set.seed(123456)
split = sample.split(data,SplitRatio = 0.7)
train_reg = subset(data, split =='TRUE')
test_reg = subset(data, split == 'FALSE')

logistic_model = glm(factor(genuine) ~ X1+X2+X3+X4+X5+X6, data = train_reg,
                     family = binomial(link = 'logit'), 
                     control = list(maxit=1000))

logistic_model
summary(logistic_model)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

predict_reg = predict(logistic_model,test_reg,type='response')
predict_reg = ifelse(predict_reg >0.5,1,0)
table(test_reg$genuine,predict_reg)
missing_classerr = mean(predict_reg!=test_reg$genuine)
print(paste('Accuracy=',1-missing_classerr))

ROCPred = prediction(predict_reg, test_reg$genuine)
ROCPer = performance(ROCPred, measure = 'tpr',x.measure = 'fpr')

auc = performance(ROCPred,measure='auc')
auc = auc@y.values[[1]]
auc

plot(ROCPer, colourise=TRUE, print.cuttoffs.at=seq(0.1,by=0.1), col=2, lwd=2, main='ROC Curve')
abline(a = 0, b = 1)
auc = round(auc, 4)
legend(.6, .4, auc, title = 'AUC', cex = 1)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

library(randomForest)

# summary
set.seed(12345)
data$genuine = as.factor(data$genuine)
split = sample.split(data,SplitRatio = 0.6)
train_rf = subset(data, split =='TRUE')
test_rf = subset(data, split == 'FALSE')
rf = randomForest(genuine~ X1+X2+X3+X4+X5+X6,data =train_rf,ntree=500)
print(rf)

1
2

importance(rf)
varImpPlot(rf)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

pred1 = predict(rf,test_rf, type='prob')
perf = prediction(pred1[,2], test_rf$genuine)
auc = performance(perf, 'auc')
auc_value = auc@y.values[[1]]

pred3 = performance(perf, 'tpr', 'fpr')

plot(pred3, main="ROC Curve for Random Forest", col=2, lwd=2)
abline(a=0, b=1, lwd=2, lty=2, col='gray')

text(0.6, 0.2, paste("AUC =", round(auc_value, 3)), col=2, cex=1.2)