キーワード言及数の統計について・・・その４

下のデータフレームはどうやって作ったの？という話。色々ゴチャゴチャやったけど、.Rhisotry見ると実際の作業はそれほど難しいものではなかったらしい。色々ゴチャゴチャ確認しながらできるのはRの大きな魅力なので、それはそれで問題なし。むしろ作業履歴をまとめてシンプル化できることがRの長所と受け止めるべし。

calcp.np <- function(v){
    dens <- density(v)
    x <- dens$x
    y <- y <- cumsum(dens$y) * (x[length(x)] - x[1])/length(x)
    p <- NULL
    for(i in v){ p <- c(p, getpercentile(x, y, i))}
    return(p)
}

getpercentile <- function(x, y, t){ # x:variable y:percentile t:point
    idx <- (which(abs(x - t) == min(abs(x - t))))[1]
    if(t == x[idx]){
        return(y[idx])
    }
    else if(t < x[idx]){
        idx <- idx - 1
    }
    return(y[idx] + (t - x[idx])*(y[idx + 1] - y[idx])/(x[idx + 1] - x[idx]))
}

df <- df[sort.list(df$total_begin, dec=TRUE),]
p1 <- calcp.np(df$total_begin)
df <- df[sort.list(df$total_3m, dec=TRUE),]
p2 <- calcp.np(df$total_3m)
df <- data.frame(df, p_3m=p2, p_begin=p1)

いちいちソートしているのは、累積分布を求めるときにソートが必要なので。関数内部でやっちゃってもよかったかなあ・・・？