差分

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

この比較画面へのリンク

次のリビジョン		前のリビジョン
mva:regression [2022/05/13 08:09] – 作成 watalumva:regression [2022/05/13 12:14] (現在) – watalu
行 1: 行 1:
 +=== 授業前に掲示した版 ===
 +
 +ビールのデータを楽しむ。
 +
 <code> <code>
 beer.1 = read.table("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_001.csv", beer.1 = read.table("http://ds.k.kyoto-u.ac.jp/e-learning_files/data_analysis_basic/jma_001.csv",
行 149: 行 153:
                  lag.max=12)                  lag.max=12)
 plot(beer.1.var) 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> </code>