Applied Statistics and Data Mining

文書の過去の版を表示しています。

ビールのデータを楽しむ。

beer.1 = read.table("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_001.csv",
                    sep=",", header=TRUE, fileEncoding="UTF-8")

# データを読み込む

beer.1 = read.csv(file("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_001.csv",
                       encoding='cp932'))
plot(beer.1$東京, beer.1$ビール, 
     main="データ:beer.1", 
     xlab="東京", 
     ylab="ビール",
     pch=20)

# データを眺める
plot(beer.1$東京)
plot(beer.1$東京, type="b")
plot(beer.1$東京, type="b", lty=2)
plot(beer.1$東京, type="b", lty=3)
plot(beer.1$東京, type="b", lty=3, pch=3)
plot(beer.1$東京, type="b", lty=3, pch=20)

plot(beer.1$ビール)
plot(beer.1$ビール, type="b")
plot(beer.1$ビール, type="b", lty=2)
plot(beer.1$ビール, type="b", lty=3)
plot(beer.1$ビール, type="b", lty=3, pch=3)
plot(beer.1$ビール, type="b", lty=3, pch=20)

hist(beer.1$東京)
hist(beer.1$東京, breaks=10)
hist(beer.1$東京, breaks=8)
hist(beer.1$東京, breaks=7)

hist(beer.1$ビール)
hist(beer.1$ビール, breaks=10)
hist(beer.1$ビール, breaks=15)
hist(beer.1$ビール, breaks=10)

boxplot(beer.1$ビール, beer.1$東京)

# データを眺める
plot(beer.1$東京, beer.1$ビール)
plot(beer.1$東京, beer.1$ビール, pch=20)
plot(beer.1$東京, beer.1$ビール, pch=20,
     xlab="東京の月平均気温", ylab="1ヶ月のビールの売上",
     main="beer.1")
plot(beer.1$東京, beer.1$ビール, pch=20,
     xlab="東京の月平均気温", ylab="1ヶ月のビールの売上",
     main="beer.1", 
     type="b")
plot(beer.1$東京, beer.1$ビール, pch=20,
     xlab="東京の月平均気温", ylab="1ヶ月のビールの売上",
     main="beer.1", 
     type="b", lty=3)

# 左上、なんかおかしい

# 月という変数を作ってみる
beer.1$月 = c(rep(c(1:12),6),c(1:4))

plot(beer.1$月, beer.1$ビール, pch=20,
     xlab="東京の月平均気温", ylab="1ヶ月のビールの売上",
     main="beer.1", 
     type="b", lty=3)

# 先ほどの左上の点が・・・
# でも今回は、この方向にモデルを頑張る話ではない
# 気温で売り上げを説明しようとすると、これは説明できない

# 元の話(回帰分析)に立ち返る

# 関数を当てはめる

lines(c(4:31), 866.168+244.7*tan((((c(4:31) - 3)/32) - 1/2) * pi/2), 
      col="red")
lines(c(4:31), 866.168+244.7*tan((((c(4:31) - 3)/32) - 1/2) * pi/2)), col="red")
lines(c(5:30), 912.8+224.6*tan((((c(5:30) - 4)/31) - 1/2) * pi/2), 
      col="red")
plot(beer.1$東京, beer.1$ビール, 
     main="データ:beer.1", 
     xlab="東京", 
     ylab="ビール",
     pch=20)
lm(ビール~東京, data=beer.1)
lines(c(0,40), c(656.93, 656.93+12.99*40))
lm(ビール~tan((((東京-4.4)/25)-1/2)*pi/2), data=beer.1)
lines(c(5:30), 877.7+188.0*tan((((c(5:30) - 4.4)/25) - 1/2) * pi/2), 
      col="red")
lm(ビール~sinh((((東京-4.4)/25)-1/2)*pi), data=beer.1)
lines(c(5:30), 878.34+86.75*sinh((((c(5:30) - 4.4)/25) - 1/2) * pi), 
      col="red")
install.packages("gam")
library("gam")
gam(ビール~東京, data=beer.1)
plot(gam(ビール~東京, data=beer.1))
summary(gam(ビール~東京, data=beer.1))
gam(ビール~s(東京), data=beer.1)
plot(gam(ビール~s(東京), data=beer.1))
beer.1.gam = gam(ビール~s(東京), data=beer.1)
plot(beer.1$東京, beer.1$ビール, 
     main="データ:beer.1", 
     xlab="東京", 
     ylab="ビール",
     pch=20)
points(beer.1$東京, predict(beer.1.gam),pch=3,col="red")

# ２乗と３乗の項を加えたい
lm(ビール~東京+東京*東京+東京*東京*東京, data=beer.1)
# 失敗

# ２乗と３乗の変数を作る
beer.1$東京2 = (beer.1$東京)^2
beer.1$東京3 = (beer.1$東京)^3
lm(ビール~東京+東京2+東京3, data=beer.1)
lines(c(0:40), c(-328.1067+257.3598*c(0:40)-16.7973*c(0:40)^2+0.3397*c(0:40)^3),
      col="blue")
# これがけっこうよさそう

# 何をしたかったんだっけ？
# 回帰分析
# 誤差を忘れてる
# 誤差が正規分布に従うモデル
# 誤差の代わりに残差を調べる

beer.1.lm3 = lm(ビール~東京+東京2+東京3, data=beer.1)
plot(beer.1.lm3)
plot(beer.1$ビール, predict(beer.1.lm3))
plot(beer.1$ビール, predict(beer.1.lm3),
     xlab="実測値", ylab="当てはめ値")
plot(fitted(beer.1.lm3), residuals(beer.1.lm3),
     xlab="当てはめ値", ylab="残差")
# 当てはめ値＋残差＝実測値
# 上の方に外れ値が残ってる

# 実は・・・

library(vars)
install.packages("vars")
library(vars)
VAR(cbind(beer.1$ビール, beer.1$東京))
VAR(cbind(beer.1$ビール, beer.1$東京), p=2)
beer.1.var = VAR(cbind(beer.1$ビール, beer.1$東京), 
                 p=2, type="const")
coef(beer.1.var)
plot(beer.1.var)
beer.1.var = VAR(cbind(beer.1$ビール, beer.1$東京), 
                 p=2, type="const",
                 lag.max=12)
plot(beer.1.var)

元の趣旨に立ち返って回帰分析を進める。

beer.2 = read.csv(file("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_002.csv",
                       encoding='cp932'))
beer.2.lm = lm(ビール~東京+京都, data=beer.2)
summary(beer.2.lm)
plot(beer.2.lm)
beer.1.lm = lm(ビール~東京, data=beer.2)
summary(beer.1.lm)
# よくなってなかった

# もう12月だけダミー変数を入れて調整してしまう
beer.3 = read.csv(file("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_003.csv",
                       encoding='cp932'))
beer.3
names(beer.3)
beer.3.lm = lm(ビール~東京+X12月, data=beer.3)
summary(beer.3.lm)
plot(beer.3.lm)
# よくなった

12月だけ異常なことを考慮したモデル

# 最後、月ダミー
beer.4 = read.csv(file("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_004.csv",
                       encoding='cp932'))
lbeer.4.lm = lm(ビール~東京+京都+東京1月後+時間+X1月+X2月+X3月+X4月+X5月+X6月+X7月+X8月+X9月+X10月+X11月+X12月+時間,data=beer.4)
summary(beer.4.lm)
plot(beer.4.lm)