授業前に掲示した版
ビールのデータを楽しむ。
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)
これをレポートにまとめてみてください。
授業中に修正した版
# データを読み込む
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")