極めて個人的なメモ

■[R]

ネットワーク分析の再勉強。

# snaパッケージ
library(sna)
el.sna <- matrix(c(
   1, 2, 1,
   2, 1, 1,
   1, 3, 1,
   3, 1, 1,
   1, 4, 1,
   4, 1, 1,
   5, 2, 1,
   5, 6, 1,
   6, 5, 1,
   6, 2, 1,
   7, 2, 1,
   8, 2, 1,
   9, 2, 1,
  10, 2, 1,
  11, 2, 1,
  12, 2, 1,
  13, 2, 1,
  14, 3, 1,
  15, 3, 1,
  16, 3, 1,
  17, 3, 1,
  18, 3, 1,
  19, 3, 1,
  20, 3, 1,
  21, 4, 1,
  22, 4, 1,
  23, 4, 1,
  24, 4, 1,
  25, 4, 1), ncol=3, byrow=TRUE)
attr(el.sna, "n") <- 25

gplot(el.sna, mode="kamadakawai", displaylabels=TRUE)
gden(el.sna)   # 密度
gtrans(el.sna) # 推移性



# igraphパッケージ
library(igraph)
el.igraph <- graph.edgelist(matrix(c(
  "A", "B",
  "B", "A",
  "A", "C",
  "C", "A",
  "A", "D",
  "D", "A",
  "E", "B",
  "F", "B",
  "E", "F",
  "F", "E",
  "G", "B",
  "H", "B",
  "I", "B",
  "J", "B",
  "K", "B",
  "L", "B",
  "M", "B",
  "N", "C",
  "O", "C",
  "P", "C",
  "Q", "C",
  "R", "C",
  "S", "C",
  "T", "C",
  "V", "D",
  "W", "D",
  "X", "D",
  "Y", "D",
  "Z", "D"), ncol=2, byrow=TRUE))
plot(el.igraph, layout=layout.kamada.kawai, vertex.label=V(el.igraph)$name)
graph.density(el.igraph) # 密度
transitivity(el.igraph)  # 推移性

ツイートする

Rのバージョンアップの際、これまで使っていたパッケージも同時にアップデートしたいんだけど、毎回その効率的なやり方を忘れてしまうのでここにメモしておく。

ネタ元は神戸大の中澤先生のメモ（ http://minato.sip21c.org/swtips/R.html ）。

1. 最新版のRをインストール
2. 旧版のRのライブラリフォルダを最新版のRのライブラリフォルダに上書き。ただし既に同名のファイルやフォルダがある場合は上書きしない。
3. 新しいバージョンのRを起動し、update.packages(checkBuilt=TRUE, ask=FALSE)を実行。

また、Rblas.dllの差し替えや、Rprofile.siteの書き換えも必要に応じて行うこと。

ツイートする

■[R]Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsj

Tokyo.SciPy#3で id:AntiBayesian さんがLTされていた「ガチャとは心の所作」が非常に興味深かったのでRでもやってみました。

（もっとも、ustがなかったので私は発表資料を読んだだけですが…）

「ガチャとは心の所作」 d:id:AntiBayesian:20120318

Rで書いたらこんな感じになりました。

ガチャのシミュレーション関数はこんな感じ。

# probが各アイテムの出現確率
# nopがガチャを回す人数（既定値は10000人）
# iter.maxが一人が回す回数の上限（既定値は10000回）
gacha.sim <- function(prob, nop=10000, iter.max=10000){
  nop.result <- numeric(nop)
  items <- length(prob)
  result <- matrix(sample(1:items, nop*iter.max, replace=TRUE, prob=prob),
                   iter.max)
  for (i in 1:nop) {
    for (j in 1:iter.max) {
      if (length(unique(result[1:j, i]))==items){
        nop.result[i] <- j
        break
        }
      }
    }
  list("result" = nop.result, "prob"=prob)
  }

二重ループで効率悪いのですが、きっと id:a_bicky さんや裏RjpWikiの中の人がもっと良い方法を示唆してくれるものと思います。

シミュレーションの実行はこんな感じで。

# ここでは４つの出現確率パターンで試してみる。
result1 <- gacha.sim(c(  1,  1,  1,  1, 1, 1))
result2 <- gacha.sim(c(100, 50, 10, 10, 3, 1))
result3 <- gacha.sim(c(  5,  5,  3,  3, 2, 1))
result4 <- gacha.sim(c( 10,  3,  2,  2, 1, 1))

シミュレーションの結果をヒストグラムにしてみましょう。

library(MASS)
windows(12, 9)
par(mfrow=c(2, 2))

x <- result1
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(30,600,
       legend=c("出現確率 [1, 1, 1, 1, 1, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

x <- result2
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(500,800,
       legend=c("出現確率 [100, 50, 10, 10, 3, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

x <- result3
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(75,600,
       legend=c("出現確率 [5, 5, 3, 3, 2, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

x <- result4
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(75,1000,
       legend=c("出現確率 [10, 3, 2, 2, 1, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

描画ももっとスマートな方法がありそう。

軸の目盛もちゃんと合わせた方がいいよなぁ。

ちなみに実行速度は

system.time(result1 <- gacha.sim(c(  1,  1,  1,  1, 1, 1)))
# ユーザ   システム       経過  
#     4.94       0.33       5.27 

system.time(result2 <- gacha.sim(c(100, 50, 10, 10, 3, 1)))
# ユーザ   システム       経過  
#    31.11       0.20      31.47 

system.time(result3 <- gacha.sim(c(  5,  5,  3,  3, 2, 1)))
# ユーザ   システム       経過  
#    5.99       0.20       6.20 

system.time(result4 <- gacha.sim(c( 10,  3,  2,  2, 1, 1)))
# ユーザ   システム       経過  
#    6.85       0.22       7.07 

ちょっと思ったこととして、我々がガキの頃１回20円で回してたガチャガチャは非復元抽出だったわけだけど、今のソーシャルゲームのガチャは復元抽出なんだよね。

そのあたりが感覚的にわかりにくくて文句言う人が出てきてるのかも。

なお、coupon collector's problemについてはこちらのびいだまブログさんの記事も参考になります。

【追記】

期待通りにビッキーさんがやってくれました！

d:id:a_bicky:20120318:1332083264

こちらに載っているgacha.sim2関数は爆速です！

一応、私の環境で試してみると、

library(rbenchmark)
p <- rep(1, 6); benchmark(gacha.sim(p), gacha.sim2(p), replications = 1)
#            test replications elapsed relative user.self sys.self user.child sys.child
# 1  gacha.sim(p)            1    6.27 2.125424      5.95     0.29         NA        NA
# 2 gacha.sim2(p)            1    2.95 1.000000      2.95     0.00         NA        NA

p <- c(100, 50, 10, 10, 3, 1); benchmark(gacha.sim(p), gacha.sim2(p), replications = 1)
#            test replications elapsed relative user.self sys.self user.child sys.child
# 1  gacha.sim(p)            1   44.72  17.0038     44.35     0.33         NA        NA
# 2 gacha.sim2(p)            1    2.63   1.0000      2.62     0.00         NA        NA

すげー。

３倍から22倍も早くなってる。

ツイートする