キーワード言及数の統計について・・・続き

で。正規分布がダメなので別な手法で行くことにする。簡単に思い浮かぶのは、正規分布じゃなくてデータ自身から累積分布関数を作ろうっていうノンパラメトリックな手法。分布関数はさっき作った密度関数から作ることにする。せっかくあるんだから使わねば。
どうやって作るか・・・というとRには累積和(?)を求めるcumsumという関数が用意されているのでそれを使いましょう。密度関数のデータはdensityの返り値のxとyという要素からいただく。

> dens <- density(df$total_3m)
> dens

Call:
        density.default(x = df$total_3m)

Data: df$total_3m (49 obs.);    Bandwidth 'bw' = 20.05

       x                y             
 Min.   :-57.14   Min.   :-6.056e-19  
 1st Qu.:200.93   1st Qu.: 1.075e-15  
 Median :459.00   Median : 7.897e-05  
 Mean   :459.00   Mean   : 9.677e-04  
 3rd Qu.:717.07   3rd Qu.: 6.676e-04  
 Max.   :975.14   Max.   : 8.559e-03

ここで問題が。densの累積和が1にならない。

> cumsum(dens$y)[length(dens$y)]
[1] 0.495457

考えてみれば当たり前で、密度関数は積分して1になることを考えると近似的にBandwidth*yの合計、y * (max(x)-min(x))/length(x)の合計を取ればよいことになります。

> x <- dens$x
> y <- cumsum(dens$y) * (x[length(x)] - x[1])/length(x)

これで累積分布関数は用意できた。次は、実際のデータ(言及数)が累積分布のどのへんにいるかという話。どのへんにいるか・・・というのを普通のプログラミング言語で作ろうと思ったら、配列を順番に見ていってx[n] < t < x[n+1]になるnを求めることになりますが、Rはwhichを使うと簡単に分かります。nを求めたら、あとは左右の点からどれだけ離れているかの配分を計って値を推定する・・・ということをしてます。

get3mpercentile <- 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
    }
    r <- (t - x[idx])/(x[idx + 1] - x[idx])
    return(y[idx] + r*(y[idx + 1] - y[idx]))
}

これで準備は完了。あとは実行するだけ。

p <- NULL
for(i in df$total_3m){ p <- c(p, get3mpercentile(x, y, i))}

今更ながらget3mpercentileを配列が取れるように作るべきだったなーと・・・。
何はともあれできましたのでplotしてみる。

一番右にいるのがらき☆すた。群抜きすぎ・・・。
こうしてみると、大体言及数200以内のアニメが全体の9割近くでダンゴになってて、それ以上の作品がいわゆる「話題作」として存在するのがわかります。言及数200以上のアニメはというと、この辺のようです。

> df[df$total_3m > 200,]
               abb total_begin total_3m
1       らき☆すた        1304      915
2     グレンラガン         585      362
3 アイドルマスター         426      305
4   ハヤテのごとく         682      300
5       電脳コイル         334      265
6   なのはStrikerS         514      236
>

もうR便利すぎ。これほど「したいこと」を簡単にやってくれるツールが他にありますか?Excelではこの軽快さはちょっと味わえないですねー。
それはそうと、結果としてはなるほどねというラインナップ。とりあえずアニオタならずとも一般常識としてこれらの作品は押さえておくべきかも。

ちょっと余談。何の水準か分からないけど、有意水準出してみましょう。

> 1-p[(1-p) < 0.05]
[1] 0.01084360 0.03129061

有意水準5%未満なのは2個。これは何かというと

> df[which((1-p) < 0.05),]
           abb total_begin total_3m
1   らき☆すた        1304      915
2 グレンラガン         585      362

この2作品。やはりというかなんというか・・・。今回の分布を、「新番組としては通常の言及数である」という帰無仮説と考えると、言及数的には「普通のアニメ」の枠を超えているのはこの2作品と言えそうです。グレンラガンも立派に話題作なんですよねー。

続いてはpercentileの比較に参ります。どうやって比較するかなーってのが難しいところ。今考えてるのが(1 - p)の比を取るというもので、これは1%を2%にするより99%を100%にするほうが難しいんじゃない?という理論に基づく方法。