# -*- mode: R -*-
##EmacsとESS使いの方に。モードラインの設定。
## 2016年10月7日(金曜日)
# -*- mode: R -*-
##EmacsとESS使いの方に、モードラインの設定。

##plot(sin)
##ディレクトリの変更は、これでもできる。
##setwd(choose.dir())


par(mfcol=c(1,2))
##表面積一定の立方体の体積は
x1 <- seq(0,2,by=0.02)
x2 <- x1
## 制約式:x1*x2+x2*x3+x3*x1=3より
## x3*(x1+x2)=3-x1*x2だから
## 体積を求める関数は、x1とx2を用いて
g <- function(x1,x2){x1*x2*(3.0-x1*x2)/(x1+x2)}
taiseki <- outer(x1,x2,g)
contour(x1,x2,taiseki)
persp(x1,x2,taiseki,theta = 30, phi = 30)

##実際に2次計画問題を解いてみる。
install.packages("e1071_1.6-7.zip", repo=NULL)
##install.packages(“e1071”, dependencies=TRUE)
##library(e1071)

##データの読み込み
source("patnew.dat")
##svmの実行
model <- svm(x, y, kernel="linear")
summary(model)
pred <- predict(model, x)

##サポートベクターに選ばれた入力データは
model$index

##kernel="linear"の場合、平面なので、ニューラルネットワークの1つのセル
##と同じ
apply(model$SV, 2, mean)

barplot(apply(-model$SV, 2, mean),
        col=rep(c("red","green","blue","yellow"),9))


#予測と本当のラベル
table(pred, y) 

##kernel関数の変更
poi <- function(k){
   model <- svm(x, y, kernel=k, coef0=1);
   summary(model);
   pred <- predict(model, x);
   table(pred, y);
   plot(cmdscale(dist(x)),
     col = as.integer(y),
     pch = c("o","+")[seq(along=y) %in% model$index + 1]);
   title(k)
}
poi("linear")
poi("polynomial") ##coef0=1, dgree=3(default)
poi("radial")

##cmdscaleは、要素間の違いに基づき、2次元平面に分布させる方法
#distは、要素間の距離行列を計算する。ユークリッド距離。
#col = as.integer(y)で、0(黒)、1(赤)で各点を書く。
#pchで、各点に現れる文字を指定する。
#model$indexは、サポートベクトルを示す
#%in%等を使って、サポートベクトルか否かの区別をしている。match(x, table)


##マージンを最大化しておくと将来役立つ。
##これが、「未学習データに対しても良い」の理屈。
maru <- matrix(
               c(13, 72,
                 11, 35,
                 37, 57,
                 23, 41,
                 26, 21,
                 40, 68,
                 43, 50,
                 68, 75,
                 89, 77),
               ncol=2 ,byrow=T)
batsu <- matrix(
                c(68, 50,
                  77, 43,
                  88, 46,
                  62, 27,
                  78, 29,
                  62, 31,
                  55, 13,
                  89, 16),
                ncol=2 ,byrow=T)

##svmのマージン最大化の説明
plot(rbind(maru,batsu),type="n")
##pch=21は色ありの●。
points(maru[,],pch=21,col=1,cex=2,lw=3,bg="orange")
points(batsu[,],pch=4,col=2,cex=2,lw=3)

col2rgb("orange")/255
plot(rbind(maru,batsu),type="n")
##21は色ありの●。
for( i in 1:1000 ){
  points(maru[,]+5*rnorm(length(as.numeric(maru))),
         col=NULL, bg=rgb(1, 0.64, 0, alpha=0.1), pch=21, cex=2)
}
points(maru[,],pch=21,col=1,cex=2,bg="orange",lw=3)

for( i in 1:1000 ){
  points(batsu[,]+5*rnorm(length(as.numeric(batsu))),
         pch=4,col=1,cex=2)
}
points(batsu[,],pch=4,col=2,cex=2,lw=2)