Hatena::ブログ(Diary)

極めて個人的なメモ Twitter

Θ・)ノシ Bob#3のメモ帳です。
1972 | 12 |
2003 | 03 | 04 | 05 | 06 | 11 | 12 |
2004 | 01 | 02 | 03 | 04 | 05 | 06 | 09 | 10 | 11 | 12 |
2005 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 09 | 10 | 11 |
2006 | 03 | 04 | 05 | 06 | 07 | 09 | 10 | 11 | 12 |
2007 | 01 | 02 | 03 | 05 | 06 | 07 | 10 | 11 | 12 |
2008 | 01 | 02 | 04 | 05 | 08 | 09 | 10 | 11 |
2009 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2010 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 |
2011 | 05 | 06 | 09 | 10 |
2012 | 02 | 03 | 05 | 06 | 08 |
2013 | 01 | 02 | 03 |
2015 | 02 |

2012-03-18

[]Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsj Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsjを含むブックマーク Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsjのブックマークコメント

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

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


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


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

f:id:bob3:20120318173714p:image:w640


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

# 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の中の人がもっと良い方法示唆してくれるものと思います。


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

# ここでは4つの出現確率パターンで試してみる。
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 

ちょっと思ったこととして、我々がガキの頃1回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

すげー。

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

r-de-rr-de-r 2012/04/13 16:20 ご指名があったことに気づきましたので,書いて見ました

r-de-rr-de-r 2012/04/15 22:42 グラフを簡単に描く方法は,明らかと思いますが...
results <- list(
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))
draw <- function(x) {
truehist(x$result, xlim = c(min(x$result), max(x$result)),
xlab = "コンプまでの回数", ylab = "人数", prob = FALSE)
legend("topright", legend = c(paste("出現確率 [", paste(x$prob,
collapse = ", "), "]"), 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 = "")
}
sapply(results, draw)

2012-03-10 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)

f:id:bob3:20120310115634p:image:w640


続いて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)

f:id:bob3:20120310151538p:image:w640