DecTheory/PlugIn.R at master · Zoncker/DecTheory · GitHub

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
## Âîññòàíîâëåíèå öåíòðà -  ìàòîæèäàíèÿ íîðìàëüíîãî ðàñïðåäåëåíèÿ
estimateMu <- function(objects)
{

  cols <- dim(objects)[2]
  mu <- matrix(NA, 1, cols) #ñîçäà¸ì âåêòîð
  for (col in 1:cols)
  {
    mu[1, col] = mean(objects[,col])
  }
  return(mu)
}
## Âîññòàíîâëåíèå êîâàðèàöèîííîé ìàòðèöû íîðìàëüíîãî ðàñïðåäåëåíèÿ
estimateCovarianceMatrix <- function(objects, mu)
{
  rows <- dim(objects)[1]
  cols <- dim(objects)[2]
  sigma <- matrix(0, cols, cols)
  for (i in 1:rows)
  {
    sigma <- sigma + (t(objects[i,] - mu) %*% # - îïåðàòîð óìíîæåíèÿ ìàòðèö, t - transpose
                        (objects[i,] - mu)) / (rows - 1)
  }
  return (sigma)
}
## Êîýôôèöèåíòû êâàäðàòè÷íîãî äèñêðèìèíàíòà
getPlugInDiskriminantCoeffs <- function(mu1, sigma1, mu2,
                                        sigma2)
{
  ## êðèâàÿ âòîðîãî ïîðÿäêà: a*x1^2 + b*x1*x2 + c*x2 + d*x1 + e*x2 + f = 0
  invSigma1 <- solve(sigma1)

  invSigma2 <- solve(sigma2)
  f <- log(abs(det(sigma1))) - log(abs(det(sigma2))) +
    mu1 %*% invSigma1 %*% t(mu1) - mu2 %*% invSigma2 %*%
    t(mu2);
  alpha <- invSigma1 - invSigma2
  a <- alpha[1, 1]
  b <- 2 * alpha[1, 2]
  c <- alpha[2, 2]
  beta <- invSigma1 %*% t(mu1) - invSigma2 %*% t(mu2)
  d <- -2 * beta[1, 1]
  e <- -2 * beta[2, 1]
  ##ðåøàåì ïðîáëåìó âîçâðàùåíèÿ èç ôóíêöèè íàáîðà ïàðàìåòðîâ ïåðåäà÷åé èõ map-îì
  return (c("x^2" = a, "xy" = b, "y^2" = c, "x" = d, "y"
            = e, "1" = f))
}
ObjectsCountOfEachClass <- 100
## Ïîäêëþ÷àåì áèáëèîòåêó MASS äëÿ ãåíåðàöèè ìíîãîìåðíîãî íîðìàëüíîãî ðàñïðåäåëåíèÿ
library(MASS)
## Ãåíåðèðóåì òåñòîâûå äàííûå
Sigma1 <- matrix(c(10, 0, 0, 1), 2, 2)
Sigma2 <- matrix(c(1, 0, 0, 5), 2, 2)
Mu1 <- c(1, 0)
Mu2 <- c(15, 0)
## Èñïîëüçóåì ôóíêöèîíàë MASS
xy1 <- mvrnorm(n=ObjectsCountOfEachClass, Mu1, Sigma1)
xy2 <- mvrnorm(n=ObjectsCountOfEachClass, Mu2, Sigma2)
## Ñîáèðàåì äâà êëàññà â îäíó âûáîðêó
xl <- rbind(cbind(xy1, 1), cbind(xy2, 2))
## Ðèñóåì îáó÷àþùóþ âûáîðêó
colors <- c(rgb(0/255, 162/255, 232/255), rgb(0/255,
                                              200/255, 0/255))
plot(xl[,1], xl[,2], pch = 21, bg = colors[xl[,3]], asp =
       1)
## Îöåíèâàíèå
objectsOfFirstClass <- xl[xl[,3] == 1, 1:2]
objectsOfSecondClass <- xl[xl[,3] == 2, 1:2]

mu1 <- estimateMu(objectsOfFirstClass)
mu2 <- estimateMu(objectsOfSecondClass)
sigma1 <- estimateCovarianceMatrix(objectsOfFirstClass,
                                   mu1)
sigma2 <- estimateCovarianceMatrix(objectsOfSecondClass,
                                   mu2)
coeffs <- getPlugInDiskriminantCoeffs(mu1, sigma1, mu2,
                                      sigma2)
## Îòðèñîâêà êðèâîé
x <- y <- seq(-10, 20, len=100)
z <- outer(x, y, function(x, y) coeffs["x^2"]*x^2 +
             coeffs["xy"]*x*y
           + coeffs["y^2"]*y^2 + coeffs["x"]*x
           + coeffs["y"]*y + coeffs["1"])
contour(x, y, z, levels=0, drawlabels=FALSE, lwd = 3, col =
          "red", add = TRUE)