Data clustering example for RBF neural network

Our data are from data06.csv

rbfdata <- read.csv("./data/data03.csv", sep = ";", header = FALSE)
colnames(rbfdata)[1] <- "x1"
colnames(rbfdata)[2] <- "x2"
colnames(rbfdata)[3] <- "d"
names(rbfdata)
## [1] "x1" "x2" "d"
ncol(rbfdata)
## [1] 3
nrow(rbfdata)
## [1] 100

Random sampling for train-set

trainset = rbfdata[sample(nrow(rbfdata), replace = F, size = 0.8 * nrow(rbfdata)), 
    ]
ncol(trainset)
## [1] 3
nrow(trainset)
## [1] 80
x = trainset$x1
xx = x
y = trainset$x2
yy = y
d = trainset$d
plot(x, y, col = (trainset$d + 1), pch = 23, cex = 1.2, xlab = "x1", ylab = "x2", 
    main = "Points from train set")
text(x + 0.01, y + 0.01, labels = as.character(1:80))

plot of chunk unnamed-chunk-2

Searching for RBF-func centers with K-means clustering (3 centroids)

dataFrame <- data.frame(xx, yy)
kmeansObj <- kmeans(dataFrame, centers = 4)
names(kmeansObj)
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"
kmeansObj$cluster
##  [1] 3 2 3 2 3 4 3 4 3 3 3 1 2 1 4 3 4 4 3 3 4 2 3 2 4 4 1 4 1 2 4 2 4 3 3
## [36] 3 3 1 1 4 3 2 3 1 2 1 4 4 4 1 4 4 1 4 2 1 1 2 2 2 4 2 1 4 3 4 3 3 1 4
## [71] 1 3 2 1 4 1 2 1 2 4

plot(xx, yy, col = d + 1, pch = 19, cex = 2)
points(kmeansObj$centers, pch = 3, cex = 3, lwd = 5)

plot of chunk unnamed-chunk-3


plot(xx, yy, col = kmeansObj$cluster, pch = 19, cex = 2)
points(kmeansObj$centers, col = kmeansObj$cluster, pch = 3, cex = 3, lwd = 5)

plot of chunk unnamed-chunk-3


kmeansObj$centers
##       xx     yy
## 1 0.3689 0.7883
## 2 0.8183 0.6600
## 3 0.6754 0.2707
## 4 0.1792 0.3337

So we have out C1 - C4 centroids for neural network:

c1 = as.vector(kmeansObj$centers[1, ])
print(c1)
## [1] 0.3689 0.7883
c2 = as.vector(kmeansObj$centers[2, ])
print(c2)
## [1] 0.8183 0.6600
c3 = as.vector(kmeansObj$centers[3, ])
print(c3)
## [1] 0.6754 0.2707
c4 = as.vector(kmeansObj$centers[4, ])
print(c4)
## [1] 0.1792 0.3337

Reprojecting our data in new with non-linear func:

phi <- function(x, c, sigm) exp(-sqrt((x[1] - c[1])^2 + (x[2] - c[2])^2)/(2 * 
    sigm^2))
h1 = x
h2 = y
h3 = x
h4 = y

for (i in 1:80) {
    X = c(x[i], y[i])
    h1[i] = phi(X, c1, 0.5)
    h2[i] = phi(X, c2, 0.5)
    h3[i] = phi(X, c3, 0.5)
    h4[i] = phi(X, c4, 0.5)

}
h1
##  [1] 0.1739 0.4453 0.3128 0.4614 0.4536 0.4828 0.2338 0.3200 0.4078 0.3731
## [11] 0.1705 0.7420 0.3645 0.8544 0.2975 0.3986 0.3148 0.3165 0.2277 0.2434
## [21] 0.4040 0.3692 0.2184 0.4270 0.2202 0.2829 0.7075 0.2086 0.8743 0.4351
## [31] 0.4096 0.2671 0.5112 0.4069 0.1944 0.3009 0.3130 0.6482 0.5511 0.4396
## [41] 0.3840 0.4955 0.5188 0.7978 0.3371 0.7144 0.3394 0.5255 0.2684 0.7349
## [51] 0.3269 0.4156 0.5368 0.4096 0.3998 0.9140 0.8697 0.4150 0.3249 0.3201
## [61] 0.4592 0.3484 0.7557 0.5112 0.1925 0.5101 0.2424 0.3892 0.6243 0.2323
## [71] 0.9538 0.3226 0.4332 0.7663 0.6011 0.5869 0.2511 0.7753 0.3160 0.2759
h2
##  [1] 0.3071 0.7607 0.3770 0.6109 0.5728 0.2857 0.3447 0.2052 0.4022 0.4335
## [11] 0.3034 0.2925 0.5790 0.3594 0.1850 0.5686 0.2312 0.2414 0.4713 0.5296
## [21] 0.2714 0.7406 0.3090 0.7638 0.1668 0.2633 0.5507 0.1649 0.3887 0.8513
## [31] 0.2812 0.5945 0.3921 0.6104 0.3418 0.3943 0.5609 0.5766 0.2280 0.2124
## [41] 0.3986 0.7924 0.4999 0.4698 0.6920 0.5358 0.2016 0.2604 0.1720 0.5262
## [51] 0.2230 0.2806 0.2186 0.2018 0.6212 0.3985 0.4394 0.5304 0.8136 0.8045
## [61] 0.2371 0.7274 0.5128 0.2336 0.2614 0.3465 0.2595 0.5341 0.2664 0.1491
## [71] 0.4063 0.4504 0.8225 0.4428 0.3129 0.2419 0.5808 0.4268 0.7887 0.1891
h3
##  [1] 0.5725 0.5428 0.7991 0.2855 0.6512 0.3798 0.7682 0.3346 0.6551 0.7666
## [11] 0.5605 0.2318 0.2547 0.3080 0.3021 0.7269 0.4062 0.4309 0.6063 0.5770
## [21] 0.4147 0.3232 0.6981 0.5554 0.3233 0.5310 0.3555 0.3291 0.2726 0.4907
## [31] 0.4301 0.5682 0.5143 0.6862 0.6357 0.8756 0.7495 0.4402 0.2275 0.2627
## [41] 0.6914 0.4253 0.5764 0.3747 0.6093 0.3381 0.3092 0.3075 0.2929 0.3893
## [51] 0.3735 0.4235 0.2139 0.2586 0.2767 0.2859 0.3445 0.2425 0.3873 0.4537
## [61] 0.3024 0.5860 0.3386 0.2638 0.5939 0.4575 0.5704 0.7619 0.1875 0.2612
## [71] 0.3009 0.9294 0.5117 0.3889 0.3521 0.1779 0.5242 0.2726 0.4623 0.3327
h4
##  [1] 0.2278 0.3107 0.4507 0.2045 0.4038 0.7733 0.3300 0.8519 0.5461 0.4688
## [11] 0.2228 0.3511 0.1651 0.4302 0.7746 0.3781 0.7930 0.7642 0.2251 0.2207
## [21] 0.8751 0.1893 0.3327 0.3050 0.5813 0.5857 0.3218 0.5437 0.3267 0.2796
## [31] 0.8428 0.2268 0.6069 0.3596 0.2461 0.4094 0.3023 0.3768 0.4576 0.6376
## [41] 0.5264 0.2771 0.4759 0.4034 0.2722 0.3111 0.8376 0.6496 0.7089 0.3720
## [51] 0.8646 0.8477 0.4309 0.6535 0.1829 0.3414 0.3946 0.1758 0.1953 0.2145
## [61] 0.7153 0.2706 0.3264 0.5760 0.3223 0.6798 0.4701 0.3939 0.2685 0.6110
## [71] 0.3600 0.3882 0.2881 0.4476 0.6170 0.2734 0.2083 0.2930 0.2155 0.7386

Weight adoptation

error = 1
# init weights
w1 = w2 = w3 = w4 = 0.01
try = 0
while (try < 100) {
    error = 0
    ww1 = ww2 = ww3 = ww4 = 0
    for (i in 1:80) {
        y_out <- w1 * h1[i] + w2 * h2[i] + w3 * h3[i] + w4 * h4[i]
        w1 = w1 - 0.01 * (y_out - d[i]) * h1[i]
        w2 = w2 - 0.01 * (y_out - d[i]) * h2[i]
        w3 = w3 - 0.01 * (y_out - d[i]) * h3[i]
        w4 = w4 - 0.01 * (y_out - d[i]) * h4[i]


    }
    # print(y_out) print(c(w1,w2,w3,w4))

    try = try + 1

}

err = 0
for (i in 1:80) {
    res = w1 * h1[i] + w2 * h2[i] + w3 * h3[i] + w4 * h4[i]
    # print(res)
    if (abs(res) < 0.5) {
        t = 0
        # print('Output is 0')
    }
    if (abs(res) > 0.5) {
        t = 1
        # print('Output is 1')
    }
    # print(d[i])
    if (d[i] != t) {
        err = err + 1
    }
}
print(paste("error on train set is ", err))
## [1] "error on train set is  4"
print(paste("level of correct classification is ", (80 - err)/80 * 100, "%"))
## [1] "level of correct classification is  95 %"
print(paste(c(w1, w2, w3, w4)))
## [1] "-0.391818566316758" "1.25878881688592"   "0.974494127950963" 
## [4] "-0.716746032453783"

Fit with least mean square

print(error)
## [1] 0
model <- lm(d ~ h1 + h2 + h3 + h4)
summary(model)
## 
## Call:
## lm(formula = d ~ h1 + h2 + h3 + h4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5948 -0.1384  0.0015  0.1138  0.5720 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.467      0.181    2.58  0.01175 *  
## h1            -0.643      0.154   -4.19  7.6e-05 ***
## h2             1.048      0.191    5.49  5.3e-07 ***
## h3             0.709      0.175    4.05  0.00012 ***
## h4            -1.014      0.174   -5.83  1.3e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.23 on 75 degrees of freedom
## Multiple R-squared:  0.801,  Adjusted R-squared:  0.791 
## F-statistic: 75.6 on 4 and 75 DF,  p-value: <2e-16
print(model$coefficients)
## (Intercept)          h1          h2          h3          h4 
##      0.4673     -0.6427      1.0475      0.7093     -1.0137
c <- as.numeric(model$coefficients)
data <- data.frame(x, y, h1, h2, h3, h4, d)
head(data, 5)
##       x     y     h1     h2     h3     h4 d
## 1 0.871 0.072 0.1739 0.3071 0.5725 0.2278 1
## 2 0.712 0.574 0.4453 0.7607 0.5428 0.3107 1
## 3 0.567 0.242 0.3128 0.3770 0.7991 0.4507 1
## 4 0.741 0.894 0.4614 0.6109 0.2855 0.2045 1
## 5 0.610 0.475 0.4536 0.5728 0.6512 0.4038 1
write.csv(data, "result/4akn.csv")
print(c(w1, w2, w3, w4))
## [1] -0.3918  1.2588  0.9745 -0.7167
err = 0
for (i in 1:80) {
    res = (c[1] + c[2] * h1[i] + c[3] * h2[i] + c[4] * h3[i] + c[5] * h4[i])
    if (abs(res) < 0.5) {
        t = 0
        # print('Output is 0')
    }
    if (abs(res) > 0.5) {
        t = 1
        # print('Output is 1')
    }
    # print(d[i])
    if (d[i] != t) {
        err = err + 1
    }
}
print(paste("error on train set is ", err))
## [1] "error on train set is  6"
print(paste("level of correct classification is ", (80 - err)/80 * 100, "%"))
## [1] "level of correct classification is  92.5 %"