=== 授業前に掲示した版 ===

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

<code>
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)
</code>

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

<code>
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)
</code>

これをレポートにまとめてみてください。

=== 授業中に修正した版 ===
<code>
# データを読み込む

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$東京)
beer.1$月数 = c(1:76)
plot(beer.1$月数, beer.1$東京)
plot(beer.1$月数, beer.1$東京, type="b")
plot(beer.1$月数, beer.1$東京, type="b", lty=2)
plot(beer.1$月数, beer.1$東京, type="b", lty=3)
plot(beer.1$月数, beer.1$東京, type="b", lty=3, pch=3)
plot(beer.1$月数, beer.1$東京, type="b", lty=3, pch=20)

plot(beer.1$ビール)
plot(beer.1$月数, beer.1$ビール, type="b")
plot(beer.1$月数, beer.1$ビール, type="b", lty=3)
plot(beer.1$月数, 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$東京, breaks=10)

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

boxplot(beer.1$ビール)
boxplot(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))
beer.1$月 = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
              1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
              1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
              1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
              1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
              1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
              1, 2, 3, 4)

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

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

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

# 関数を当てはめる
plot(beer.1$東京, beer.1$ビール, 
     main="データ:beer.1", 
     xlab="東京", 
     ylab="ビール",
     pch=20)
lm(ビール~東京, data=beer.1)
summary(lm(ビール~東京, data=beer.1))
lines(c(0,40), c(656.931, 656.93+12.992*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="blue")
legend("bottomright",
       legend=c("直線",  "tan", "sinh"), 
       col=c("black", "red", "blue"),
       lty=1)

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="green")
legend("bottomright",
       legend=c("直線",  "tan", "sinh", "gam"), 
       col=c("black", "red", "blue", "green"),
       lty=1)

# ２乗と３乗の項を加えたい
lm(ビール~東京+東京*東京+東京*東京*東京, data=beer.1)
lm(ビール~東京, data=beer.1)
# 失敗
# Rのモデル式の掛け算記号*は、掛け算の意味ではない

# ２乗と３乗の変数を作る
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="orange")
# これがけっこうよさそう
legend("bottomright",
       legend=c("直線",  "tan", "sinh", "gam", "3次関数"), 
       col=c("black", "red", "blue", "green", "orange"),
       lty=1)


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

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月だけ異常なことを考慮したモデル

plot(beer.1$東京, beer.1$ビール, 
     main="データ:beer.1", 
     xlab="東京", 
     ylab="ビール",
     pch=20)
points(beer.3$東京, predict(beer.3.lm),pch=3,col="green")


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