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)