ggplot2の練習中
ggplot2は記法が独特なので敬遠してきたのですが、最近ggplot2を使うと良さそうな事案に遭遇したので勉強中です。
library(ggplot2) p1 <- ggplot(iris, aes(Petal.Length, fill=Species)) p1 + geom_histogram() + opts(title="積み重ねヒストグラム") p1 + geom_histogram() + facet_grid(Species~.) + opts(title="群別ヒストグラム") p1 + geom_freqpoly(aes(color=Species)) + opts(title="度数分布") p1 + geom_density(alpha=0.7) + opts(title="密度推定") ggplot(iris, aes(x=Species, y=Petal.Length)) + geom_boxplot() + geom_jitter(aes(color=Species)) + opts(title="箱ひげ図") ggplot(iris, aes(x=Petal.Length, y=Petal.Width)) + geom_density2d(aes(color=Species)) + geom_point(aes(color=Species, alpha=0.7)) + opts(title="二次元密度推定") + scale_color_brewer(palette = "Dark2")
■
ネットワーク分析の再勉強。
# 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のバージョンアップ時のパッケージのアップデート方法
Rのバージョンアップの際、これまで使っていたパッケージも同時にアップデートしたいんだけど、毎回その効率的なやり方を忘れてしまうのでここにメモしておく。
ネタ元は神戸大の中澤先生のメモ（ http://minato.sip21c.org/swtips/R.html ）。
- 最新版のRをインストール
- 旧版のRのライブラリフォルダを最新版のRのライブラリフォルダに上書き。ただし既に同名のファイルやフォルダがある場合は上書きしない。
- 新しいバージョンのRを起動し、update.packages(checkBuilt=TRUE, ask=FALSE)を実行。
また、Rblas.dllの差し替えや、Rprofile.siteの書き換えも必要に応じて行うこと。
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倍も早くなってる。
RでPSM分析 #TokyoR
まずは直接測定法
# テスト用データ cheap <- c(rep(2, 5), rep(4, 20), rep(6, 25), rep(8, 30), rep(10, 15), rep(12, 5)) expensive <- c(rep(12, 5), rep(14, 5), rep(16, 15), rep(18, 25), rep(20, 30), rep(22, 15), rep(24, 5)) # 度数の算出、テーブルをデータフレームに変換しマージする A.df <- data.frame(table(cheap)) B.df <- data.frame(table(expensive)) colnames(A.df) <- c("PRICE", "A.freq") colnames(B.df) <- c("PRICE", "B.freq") AB.df <- merge(A.df, B.df, all=TRUE) AB.df[is.na(AB.df)] <- 0 # 価格順にソート AB.df[, 1] <- type.convert(as.character(AB.df[, 1])) AB.df <- AB.df[order(AB.df[, 1]), ] # 累積和とその反転、Aは安い方が100％、Bは高い方が100％ AB.df[, 2] <- rev(cumsum(rev(AB.df[, 2]))) AB.df[, 3] <- cumsum(AB.df[, 3]) N <- length(cheap) for (i in 2:3) { AB.df[, i+2] <- N - AB.df[,i] } AB.df[, 2:5] <- AB.df[, 2:5]/N*100 ###交点######################################################################## # 2直線の交点座標を求める関数 - 裏 RjpWiki # <http://blog.goo.ne.jp/r-de-r/e/e81316c26c521b31e9d47526f9bd5861> # (x1,y1)と(x2,y2)を通る直線と，(x3,y3)と(x4,y4)を通る直線の交点の座標を求める # xだけ算出するように改造 ############################################################################### prog <- function(x1, y1, x2, y2, x3, y3, x4, y4) { a1 <- (y2-y1)/(x2-x1) a3 <- (y4-y3)/(x4-x3) x <- (a1*x1-y1-a3*x3+y3)/(a1-a3) y <- (y2-y1)/(x2-x1)*(x-x1)+y1 return(c(x, y)) } # 最低価格 GAP <- AB.df[,2] - AB.df[,4] if (max(GAP[GAP <= 0])==0){ MIN <- c(AB.df[GAP==0,1],50) } else { X <- AB.df[, c(1,2,4)][GAP == max(GAP[GAP<0]), ] X <- X[X$PRICE==min(X$PRICE),] Y <- AB.df[, c(1,2,4)][GAP == min(GAP[GAP>0]), ] Y <- Y[Y$PRICE==max(Y$PRICE),] MIN <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } MIN.p <- round(as.numeric(MIN[1]), 0) # 最高価格 GAP <- AB.df[,3] - AB.df[,5] if (max(GAP[GAP <= 0])==0){ MAX <- c(AB.df[GAP==0,1],50) } else{ X <- AB.df[, c(1,3,5)][GAP == max(GAP[GAP<0]), ] X <- X[X$PRICE==max(X$PRICE),] Y <- AB.df[, c(1,3,5)][GAP == min(GAP[GAP>0]), ] Y <- Y[Y$PRICE==min(Y$PRICE),] MAX <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } MAX.p <- round(as.numeric(MAX[1]), 0) # 最適価格 GAP <- AB.df[,2] - AB.df[,3] if (max(GAP[GAP <= 0])==0){ OPT <- c(AB.df[GAP==0,1:2]) } else{ X <- AB.df[, c(1,2,3)][GAP == max(GAP[GAP < 0]), ] X <- X[X$PRICE==min(X$PRICE),] Y <- AB.df[, c(1,2,3)][GAP == min(GAP[GAP > 0]), ] Y <- Y[Y$PRICE==max(Y$PRICE),] OPT <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } OPT.p <- round(as.numeric(OPT[1]), 0) # 価格幅 RANGE <- MAX.p - MIN.p # 描画 windows(16,9) matplot(AB.df[,1], AB.df[,2:5], type="l", lwd=2, lty=1:2, col=c("#0041FF", "#FF3300"), main="直接測定法", xlab = "（ドル）", ylab = "％", las=2, xlim=c(min(cheap),max(expensive))) abline(h=50, col = "#728490", lty="dashed") abline(v=seq(min(AB.df[,1]), max(AB.df[,1]), 1), lty="dotted", col = "#E0E0E0") # 指標の点 points(MAX[1], MAX[2], pch=21, lwd=2, bg="#FF9999", cex=1.25) points(OPT[1], OPT[2], pch=22, lwd=2, bg="#FFFF00", cex=1.25) points(MIN[1], MIN[2], pch=23, lwd=2, bg="#66CCFF", cex=1.25) # 累積比率の凡例 matplotのデフォは（lty = 1:5, col = 1:6） legend(locator(1), lty=1:3, col=c("#0041FF", "#FF3300"), lwd=1.5, bg="white", cex=0.75, legend=c("安い", "高い")) # 指標の表示 TEMP <- c(N, MAX.p, OPT.p, MIN.p, RANGE) TEMP <- format(TEMP, width=max(nchar(TEMP))) LEGEND <- c(paste(" n=", TEMP[1], sep=""), paste("上限価格: ", TEMP[2], "ドル", sep=""), paste("最適価格: ", TEMP[3], "ドル", sep=""), paste("下限価格: ", TEMP[4], "ドル", sep=""), paste(" RANGE: ", TEMP[5], "ドル", sep="")) legend(locator(1), legend=LEGEND, pch=c(NA,21,22,23,NA), pt.bg=c(NA, "#FF9999", "#FFFF00", "#66CCFF", NA), bg="white", cex=0.75)
続いてPSM。バグが残ってないことを祈るのみ…
######################################## # PSM分析 van Westendorpの方式 ######################################### cheap <- c(rep(8,5), rep(10,20), rep(12,25),rep(14,30),rep(16,15), rep(18,5)) expensive <- c(rep(6, 5), rep(8, 5), rep(10, 15), rep(12, 25), rep(14, 30), rep(16, 15), rep(18, 5)) too_expensive <- c(rep(12, 5), rep(14, 5), rep(16, 15), rep(18, 25), rep(20, 30), rep(22, 15), rep(24, 5)) too_cheap <- c(rep(2, 5), rep(4, 20), rep(6, 25), rep(8, 30), rep(10, 15), rep(12, 5)) # A いくら以下だったら「安い」cheap # B いくら以上だったら「高い」expensive # C いくら以上だったら「高すぎる」too expensive # D いくら以下だったら「安すぎる」too cheap # 度数の算出、テーブルをデータフレームに変換 A.df <- data.frame(table(cheap)) B.df <- data.frame(table(expensive)) C.df <- data.frame(table(too_expensive)) D.df <- data.frame(table(too_cheap)) # 行に名前を付ける colnames(A.df) <- c("PRICE", "A.freq") colnames(B.df) <- c("PRICE", "B.freq") colnames(C.df) <- c("PRICE", "C.freq") colnames(D.df) <- c("PRICE", "D.freq") # マージする AB.df <- merge(A.df, B.df, all=TRUE) CD.df <- merge(C.df, D.df, all=TRUE) ABCD.df <- merge(AB.df, CD.df, all=TRUE) # NAを0に置き換え(20101219) ABCD.df[is.na(ABCD.df)] <- 0 # 価格順にソート ABCD.df[,1] <- type.convert(as.character(ABCD.df[,1])) ABCD.df <- ABCD.df[order(ABCD.df[, 1]), ] # 累積和とその反転 #Aは安い方が100％ ABCD.df[,2] <- rev(cumsum(rev(ABCD.df[,2]))) #Bは高い方が100％ ABCD.df[,3] <- cumsum(ABCD.df[,3]) #Cは高い方が100％ ABCD.df[,4] <- cumsum(ABCD.df[,4]) #Dは安い方が100％ ABCD.df[,5] <- rev(cumsum(rev(ABCD.df[,5]))) #それぞれの反転 N <- length(cheap) for (i in 2:5) { ABCD.df[, i+4] <- N - ABCD.df[,i] } ABCD.df[, 2:9] <- ABCD.df[, 2:9]/N*100 # A1〜D1がそのままの累積比率（％）、A2〜D2が反転させた数値 colnames(ABCD.df) <- c("PRICE", "A1", "B1", "C1", "D1", "A2", "B2", "C2", "D2") #A1[,2]:Cheap ☆ #B1[,3]:Expensive ☆ #C1[,4]:Too Expensive to Buy ★☆ #D1[,5]:Too Cheap to Buy ★☆ #A2[,6]:Not Feel Cheap ★ #B2[,7]:Not Feel Expensive ★ #C2[,8] #D2[,9] ###交点######################################################################## # 2直線の交点座標を求める関数 - 裏 RjpWiki # <http://blog.goo.ne.jp/r-de-r/e/e81316c26c521b31e9d47526f9bd5861> # (x1,y1)と(x2,y2)を通る直線と，(x3,y3)と(x4,y4)を通る直線の交点の座標を求める # xだけ算出するように改造 ############################################################################### prog <- function(x1, y1, x2, y2, x3, y3, x4, y4) { a1 <- (y2-y1)/(x2-x1) a3 <- (y4-y3)/(x4-x3) x <- (a1*x1-y1-a3*x3+y3)/(a1-a3) y <- (y2-y1)/(x2-x1)*(x-x1)+y1 return(c(x, y)) } ### 交点の算出 # 「安くない」「高くない」式 # 下限価格(MCP, Marginal Cheap Point) # D1[,5]:Too Cheap to Buy ★☆ # A2[,6]:Not Feel Cheap ★ GAP <- ABCD.df[,5] - ABCD.df[,6] if (max(GAP[GAP <= 0])==0){ MCP.P <- c(ABCD.df[GAP==0,c(1,5)]) } else { X <- ABCD.df[, c(1,5,6)][GAP == max(GAP[GAP<0]), ] X <- X[X$PRICE==min(X$PRICE),] Y <- ABCD.df[, c(1,5,6)][GAP == min(GAP[GAP>0]), ] Y <- Y[Y$PRICE==max(Y$PRICE),] MCP.P <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } MCP <- round(as.numeric(MCP.P[1]), 0) #上限価格(MEP, Marginal Expensive Point) # C1[,4]:Too Expensive to Buy ★☆ # B2[,7]:Not Feel Expensive ★ GAP <- ABCD.df[,4] - ABCD.df[,7] if (max(GAP[GAP <= 0])==0){ MEP.P <- c(ABCD.df[GAP==0,c(1,4)]) } else { X <- ABCD.df[, c(1,4,7)][GAP == max(GAP[GAP<0]), ] X <- X[X$PRICE==min(X$PRICE),] Y <- ABCD.df[, c(1,4,7)][GAP == min(GAP[GAP>0]), ] Y <- Y[Y$PRICE==max(Y$PRICE),] MEP.P <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } MEP <- round(as.numeric(MEP.P[1]), 0) #最適価格(OPP, Optimum Pricing Point) # C1[,4]:Too Expensive to Buy ★☆ # D1[,5]:Too Cheap to Buy ★☆ GAP <- ABCD.df[,4] - ABCD.df[,5] if (max(GAP[GAP <= 0])==0){ OPP.P <- c(ABCD.df[GAP==0,c(1,4)]) } else { X <- ABCD.df[, c(1,4,5)][GAP == max(GAP[GAP<0]), ] X <- X[X$PRICE==min(X$PRICE),] Y <- ABCD.df[, c(1,4,5)][GAP == min(GAP[GAP>0]), ] Y <- Y[Y$PRICE==max(Y$PRICE),] OPP.P <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } OPP <- round(as.numeric(OPP.P[1]), 0) #妥協価格(IDP, Indifferent Point) # A2[,6]:Not Feel Cheap ★ # B2[,7]:Not Feel Expensive ★ GAP <- ABCD.df[,6] - ABCD.df[,7] if (max(GAP[GAP <= 0])==0){ IDP.P <- c(ABCD.df[GAP==0,c(1,6)]) } else { X <- ABCD.df[, c(1,6,7)][GAP == max(GAP[GAP<0]), ] X <- X[X$PRICE==min(X$PRICE),] Y <- ABCD.df[, c(1,6,7)][GAP == min(GAP[GAP>0]), ] Y <- Y[Y$PRICE==max(Y$PRICE),] IDP.P <- prog(X[1], X[2], Y[1], Y[2], X[1], X[3], Y[1], Y[3]) } IDP <- round(as.numeric(IDP.P[1]), 0) RANGE <- MEP - MCP ############################################################################### #作図 ############################################################################### windows(16,9) # 「安くない」「高くない」式 matplot(ABCD.df[,1], ABCD.df[,4:7], type="l", pch=c(21,22,23,24), lwd=2, col=c("#FF3300", "#990066", "#339966", "#0041FF"), main="van Westendorp \n Price Sensitivity Meter", xlab = "（ドル）", ylab = "％", las=2, xlim=c(min(too_cheap),max(too_expensive))) abline(h=50, col = "#728490", lty="dashed") abline(v=seq(min(ABCD.df[,1]), max(ABCD.df[,1]), 1), lty="dotted", col = "#E0E0E0") # 指標の点 points(MEP.P[1], MEP.P[2], pch=21, lwd=2, bg="#FF9999", cex=1.25) points(IDP.P[1], IDP.P[2], pch=22, lwd=2, bg="#FFFF00", cex=1.25) points(OPP.P[1], OPP.P[2], pch=23, lwd=2, bg="#66CCFF", cex=1.25) points(MCP.P[1], MCP.P[2], pch=24, lwd=2, bg="#FF9900", cex=1.25) # 累積比率の凡例 matplotのデフォは（lty = 1:5, col = 1:6） legend(locator(1), lty=1:4, col=c("#FF3300", "#990066", "#339966", "#0041FF"), lwd=1.5, bg="white", cex=0.75, legend=c("高すぎる", "安すぎる", "安くない", "高くない")) ########################### # 指標の表示 TEMP <- c(N, MEP, IDP, OPP, MCP, RANGE) TEMP <- format(TEMP, width=max(nchar(TEMP))) LEGEND <- c(paste(" n=", TEMP[1], ".", sep=""), paste("上限価格: ", TEMP[2], "ドル", sep=""), paste("妥協価格: ", TEMP[3], "ドル", sep=""), paste("最適価格: ", TEMP[4], "ドル", sep=""), paste("下限価格: ", TEMP[5], "ドル", sep=""), paste(" RANGE: ", TEMP[6], "ドル", sep="")) legend(locator(1), legend=LEGEND, pch=c(NA,21,22,23,24,NA), pt.bg=c(NA, "#FF9999", "#FFFF00", "#66CCFF","#FF9900", NA), bg="white", cex=0.75)
Rでコレスポンデンス分析
ちょっとリクエストがあったのでRでコレスポンデンス分析（対応分析）の例。
メジャーなのはMASSパッケージに入っているcorresp関数を使うもので『S-PLUSによる統計解析』『Rによる統計解析』『R言語逆引きハンドブック』などでも紹介されています。金先生のページ(http://mjin.doshisha.ac.jp/R/26/26.html)でも解説されてますね。
でも個人的におすすめなのはcaパッケージのca関数を使う方法。英語の解説はhttp://www.statmethods.net/advstats/ca.html、日本語ではhttp://eau.uijin.com/advstats/ca.html。
- 作者: W.N.ヴェナブルズ,B.D.リプリー,W.N. Venables,B.D. Ripley,伊藤幹夫,戸瀬信之,大津泰介,中東雅樹
- 出版社/メーカー: シュプリンガー・フェアラーク東京
- 発売日: 2001/07
- メディア: 単行本
- クリック: 9回
- この商品を含むブログ (7件) を見る
- 作者: 青木繁伸
- 出版社/メーカー: オーム社
- 発売日: 2009/04/01
- メディア: 単行本
- 購入: 10人 クリック: 123回
- この商品を含むブログ (34件) を見る
- 作者: 石田基広
- 出版社/メーカー: シーアンドアール研究所
- 発売日: 2012/01/26
- メディア: 単行本（ソフトカバー）
- 購入: 4人 クリック: 61回
- この商品を含むブログ (11件) を見る
両方、簡単に実行方法を載せておきます。
まずは分析する分割表の入力。
# 分析するクロス集計表はこちらから拝借 # http://www.geocities.co.jp/WallStreet/7166/correspo/corr_02.html DAT <- matrix(c(41,23, 5, 5, 2, 15, 4, 0, 2, 1, 7,11, 2, 0, 0, 3, 3, 0, 0, 0), 4, byrow=TRUE) rownames(DAT) <- c("会社員","自営業","学生","その他") colnames(DAT) <- c("マンション","ハイツ・コーポ","アパート","テラスハウス","一戸建て")
MASSパッケージのcorresp関数を使った例。
library(MASS) DAT.ca1 <- corresp(DAT, nf=2) DAT.ca1 # 結果の表示 biplot(DAT.ca1, main="corresp()を使った例") # 描画
次にcaパッケージのca関数を使った例。
library(ca) DAT.ca2 <- ca(DAT) DAT.ca2 # 結果の表示 plot(DAT.ca2, main="ca()を使った例１") plot(DAT.ca2, mass = TRUE, contrib = "absolute", map ="rowgreen", arrows = c(FALSE, TRUE), main="ca()を使った例2")
まあ、大体こんな感じです。
