差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
mva:regression [2022/05/13 08:21] – watalumva:regression [2022/05/13 12:14] (現在) – watalu
行 1: 行 1:
 +=== 授業前に掲示した版 ===
 +
 ビールのデータを楽しむ。 ビールのデータを楽しむ。
  
行 186: 行 188:
  
 これをレポートにまとめてみてください。 これをレポートにまとめてみてください。
 +
 +=== 授業中に修正した版 ===
 +<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>