Rでカーネル密度推定-識別モデル

勉強した内容の備忘録
**今回の記事では、画像によって使用している乱数データが異なるため、すべてのデータで図形が一致していない場合があります**
****
データの作成
対応のあるデータと仮定する

```R
m10<-rnorm(1000,10,1)
m25<-rnorm(1000,25,3)
```

どんなデータかplotしてみる。

```R
#http://cse.naro.affrc.go.jp/takezawa/r-tips/r/55.html

#図を描く部分をlayoutで指定してやる
#2×2のプロットエリアをつくる。
#左上に図1,左下に図2,右下に図3が入力される
layout(matrix(c(1,0,2,3),2,2,byrow=TRUE),c(2,1),c(1,2),TRUE)

m_data=data.frame(m10=m10,m25=m25)

#図1
par(mar=c(0,3,1,1))#図の間の隙間を狭くしてやる
hist_m10_data <- hist(m10,right=FALSE,plot=FALSE)
barplot(hist_m10_data$counts,space=0,col=0)

#図2
par(mar=c(3,3,1,1))
plot(m_data,col=densCols(m10,m25),pch=5,cex=1)

library(ellipse)#楕円を3σ区間まで書き込む
var_m10_m25 <- var(m_data)
mean_m10 <- mean(m_data$m10)
mean_m25 <- mean(m_data$m25)
lines(ellipse(var_m10_m25, level=0.6827, centre=c(mean_m10,mean_m25), add=T),lwd=1,lty=1,col=2)
lines(ellipse(var_m10_m25, level=0.9545, centre=c(mean_m10,mean_m25), add=T),lwd=1,lty=1,col=2)
lines(ellipse(var_m10_m25, level=0.9973, centre=c(mean_m10,mean_m25), add=T),lwd=1,lty=1,col=2)

#図3
par(mar=c(3,0,1,1))
hist_m25_data <- hist(m25,right=FALSE,plot=FALSE)
barplot(hist_m25_data$counts,space=0,col=0,horiz=T)
```

![0011.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/db92525c-90ca-456f-41e2-bf1b7eeb72e6.png)



plotとhistでデータの密度をなんとなく確認できた。


#密度推定

```
#https://qiita.com/qooa/items/93a5410be54af9031a27
#http://cse.naro.affrc.go.jp/takezawa/r-tips/r/61.html
```

MASSから二次元カーネル密度推定を行うkde2dを呼び出す

```R
library("MASS")
dens <- kde2d(m10,m25,n=50)
#nで次元を指定
#dim(dens$z)
#[1] 50,50
```

この計算された密度を確認してみる。

```R
image(dens)
```
![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0d2d9b51-bcdc-85c9-3622-50516e84b3fa.png)

imgで表示すると、画像で密度の高い部分が確認できる。

もう少しわかりやすい図示

```R
library(epade)
par(mfrow=c(1,1))
par(mar=c(5,1,5,1))
bar3d.ade(dens$z,xw=5,zw=0.1,axes=FALSE,col="pink",wall=3,alpha=1)
```

![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/8198e81e-302c-7eb6-286b-8110cfdbd73f.png)

3Dで確認できた。
等高線でも確認してみる

```R
#データを切りのいい段階で分けるpretty
#密度の区切りを10個の区切りで切って、色を変えることで密度の高い部分をわかりやすくする
lev <- pretty(range(dens$z),10)
filled.contour(dens,levels=lev,xlim=c(5,15),ylim=c(15,35))
filled.contour(dens, levels=lev,col=gray((length(lev):1)/length(lev)),xlim=c(5,15),ylim=c(15,35))
```
![05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/86b901bc-8846-402b-8200-a4deb96857de.png)

![06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/872a1074-6213-215d-02a3-1b80bf8e6d0f.png)

もっと詳細に等高線を使って確認する。
95%について書いてみる

```R
library(spam)
library(maps)
library(fields)

#引数に指定したもののすべての組み合わせをつくるexpand.grid
#http://cse.naro.affrc.go.jp/takezawa/r-tips/r/46.html
grid  <- expand.grid(dens$x,dens$y)
grid2 <- cbind(grid, data.frame(as.vector(dens$z)))

#https://www.rdocumentation.org/packages/fields/versions/9.6/topics/Tps
#http://www.okadajp.org/RWiki/?%E7%9F%A5%E3%81%A3%E3%81%A6%E3%81%84%E3%82%8B%E3%81%A8%E3%81%84%E3%81%A4%E3%81%8B%E5%BD%B9%E3%81%AB%E7%AB%8B%E3%81%A4%28%3F%29%E9%96%A2%E6%95%B0%E9%81%94
#http://www.math.keio.ac.jp/~kei/GDS/2nd/spline.html

#thin plate spline法で計算

model  <- Tps(grid2[,1:2],grid2[,3])

#作成したmodelをpredictに当てはめて予想値を入れる
m_data<-data.frame(m10=m10,m25=m25)
m_data$dens <- predict(model, m_data[,1:2])

#http://www.okadajp.org/RWiki/?%E3%82%B0%E3%83%A9%E3%83%95%E3%82%A3%E3%83%83%E3%82%AF%E3%82%B9%E5%8F%82%E8%80%83%E5%AE%9F%E4%BE%8B%E9%9B%86%EF%BC%9A%E3%82%A4%E3%83%A1%E3%83%BC%E3%82%B8%E5%9B%B3

#データ中で0.005を満たさないものをplotする
plot(m_data[m_data$dens<0.005,1:2],xlim=c(5,15),ylim=c(15,35),pch=16)
contour(dens,add=T)
```
![101.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e428db61-2136-0a79-6282-b5bce8925f1d.png)


```
nrow(m_data[m_data$dens<0.005,1:2])
```

で確認してみると、1000個中85個のデータが外れていた。
測定データなどでは、これらを外れ値と扱うのかもしれない。

ggplotを使って95%を描く場合

```R

library(ggplot2)

#https://stackoverflow.com/questions/23437000/how-to-plot-a-contour-line-showing-where-95-of-values-fall-within-in-r-and-in

getLevel <- function(x,y,prob=0.95) {
    kk <- dens
    dx <- diff(dens$x[1:2])
    dy <- diff(dens$y[1:2])
    sz <- sort(dens$z)
    c1 <- cumsum(sz) * dx * dy
    approx(c1, sz, xout = 1 - prob)$y
}

AA<-data.frame(m10=m10,m25=m25)
L95 <- getLevel(AA$m10,AA$m25)

ggplot(AA,aes(m10,m25)) +
   stat_density2d(geom="tile", aes(fill = ..density..), contour = FALSE)+
   stat_density2d(colour="red",breaks=L95)
```
![100.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4ad0e598-7518-bdd3-2f96-7f74079da910.png)


95%を外れていたものを赤色にしてplot

```R
par(mfrow=c(2,1))
plot(m_data$m10,col=c("blue","red")[unclass(m_data$result+1)],pch=20)
plot(m_data$m25,col=c("blue","red")[unclass(m_data$result+1)],pch=20)
```

![200.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/68d9f421-56bb-a499-59c7-b9f7f99067a3.png)