廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Rcpp: カーネル密度推定のバンド幅を一個抜き交差検証法で決める

今日の川柳

バンド幅 h をいろいろ変えて一個抜き交差検証法で評価した対数尤度が結構なめらかな形になったのでブレント法で最適なバンド幅を選んでみた。

対象としたデータはこれ。

f:id:abrahamcow:20180320011441p:plain

h をいろいろ変えて一個抜き交差検証法で評価した対数尤度のプロットはこちら。

f:id:abrahamcow:20180320011532p:plain

選ばれた最適なバンド幅で推定した密度関数がこんな感じ。

f:id:abrahamcow:20180320011656p:plain

Rcpp のコード:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double kernel_logistic(double y, NumericVector x, double bw){
  return mean(1/(exp((y-x)/bw)+exp((x-y)/bw)+2))/bw;
}

Rcpp::NumericVector omit(NumericVector x, int i) {
  x.erase(i);
  return x;
}

// [[Rcpp::export]]
double looll(double bw, NumericVector x) {
  int N = x.length();
  double ll=0;
  for(int i=0;i<N;i++){
    ll =ll+ log(kernel_logistic(x[i],omit(x,i),bw));
  }
  return ll;
}

R のコード:

library(Rcpp)
sourceCpp("kernel_logis.cpp")
x <- faithful$waiting
hist(x,breaks = "FD")
h <- seq(0,2,by=0.01)
sap<-sapply(h,function(h)looll(bw=h,x=x))
plot(h,sap,ylab = "log-likelihood")
opt_looll <-optim(0.2,looll,x=x,control=list(fnscale=-1),method = "Brent",
                  lower = 0, upper=2)
kernel_logistic0 <- function(x,bw){
  function(y){sapply(y,function(y){kernel_logistic(y,x,bw)})}
}
kernel_logistic_x <- kernel_logistic0(x,opt_looll$par)

hist(x,freq = FALSE,breaks = "FD")
curve(kernel_logistic_x(x),add = TRUE,col="royalblue",lwd=2)

参考:カーネル (統計学) - Wikipedia