メモ：地震の分析

天気.jpから取得した地震データを分析したプログラムのメモです。コピペで動作するはずです。

library(plyr)

eq  <- read.csv("http://www.ianalysisllc.com/app/download/4590900767/eq_en.csv?t=1313125288", as.is=T)
eq2 <- ddply(eq, .(date2, area2), summarise, count=sum(!is.na(date1)))

eq_glm1 <- glm(count ~ area2, data = eq2, family = poisson())
eq_glm2 <- glm(count ~ area2, data = eq2, family = quasipoisson())
summary(eq_glm1)
summary(eq_glm2)

eq2$date3 <- as.Date(as.character(eq2$date2), "%Y-%m-%d") - as.Date("20110101", "%Y%m%d")
eq2       <- eq2[order(eq2$date3), ]
head(eq2)


plot(eq2$date3[eq2$area2=="Touhoku"], eq2$count[eq2$area2=="Touhoku"], type="l")
plot(eq2$date3[eq2$area2=="Kanto"], eq2$count[eq2$area2=="Kanto"], type="l")


eq_glm3 <- glm(count ~ area2 + date3, data = eq2, family = quasipoisson())
summary(eq_glm3)

exp(1.02 + 1.25)
exp(1.02 - 0.93)
exp(1.02 + 1.25 + 100*0.00159)

偶然にもログが残っていたので、いろいろ試行錯誤しながらエラーを出しまくってるプログラムも残しておきます。

library(plyr)
eq  <- read.csv("http://www.ianalysisllc.com/app/download/4590900767/eq_en.csv?t=1313125288", as.is=T)
ddply(eq, .(date1), summarise, count=sum(!is.na(date2)))
ddply(eq, .(date1, area), summarise, count=sum(!is.na(date2)))
str(eq)
ddply(eq, .(date1, area2), summarise, count=sum(!is.na(date2)))
peq_glm <- glm(count ~ area2, data = eq2, family = quasipoisson())
eq2 <- ddply(eq, .(date1, area2), summarise, count=sum(!is.na(date2)))
peq_glm <- glm(count ~ area2, data = eq2, family = quasipoisson())
eq_glm <- glm(count ~ area2, data = eq2, family = quasipoisson())
summary(eq_glm)
str(eq2)
eq2$date3 <- eq2$date2 - 20110101
str(eq2)
eq2$date3 <- eq2$date1 - 20110101
str(eq2)
eq_glm <- glm(count ~ area2, data = eq2[226:755, ], family = quasipoisson())
eq_glm <- glm(count ~ area2 + date3, data = eq2[226:755, ], family = quasipoisson())
summary(eq_glm)
plot(eq2$date3, eq2$count)
plot(eq2$date3[226:755], eq2$count[226:755])
eq2$date3 <- as.Date(as.character(eq2$date1), "%Y%m#d") - as.Date("20110101", "%Y%m#d")
head(eq2)
tail(eq2)
as.Date(as.character(eq2$date1), "%Y%m#d")
as.Date(as.character(eq2$date1), "%Y%m%d")
eq2$date3 <- as.Date(as.character(eq2$date1), "%Y%m%d") - as.Date("20110101", "%Y%m#d")
tail(eq2)
eq2$date3 <- as.Date(as.character(eq2$date1), "%Y%m%d") - as.Date("20110101", "%Y%m%d")
tail(eq2)
plot(eq2$date3[226:755], eq2$count[226:755])
eq2[300:400,]
eq2
eq2 <- ddply(eq, .(date2, area2), summarise, count=sum(!is.na(date1)))
head(eq2)
eq2$date3 <- as.Date(as.character(eq2$date1), "%Y-%m-%d") - as.Date("20110101", "%Y%m%d")
eq2$date3 <- as.Date(as.character(eq2$date2), "%Y-%m-%d") - as.Date("20110101", "%Y%m%d")
head(eq2)
plot(eq2$date3, eq2$count)
plot(eq2$date3, eq2$count, type="l")
plot(eq2$date3[eq2$area2=="Touhoku"], eq2$count[eq2$area2=="Touhoku"], type="l")
eq2 <- eq2[order(eq2$date3), ]
plot(eq2$date3[eq2$area2=="Touhoku"], eq2$count[eq2$area2=="Touhoku"], type="l")
eq_glm <- glm(count ~ area2 + date3, data = eq2, family = quasipoisson())
summary(eq_glm)
(a <- exp(4.529-0.03883*20))
exp(1.02+1.25)
exp(1.02-0.93)
exp(1.02+1.25 + 100*0.00159)
eq_glm2 <- glm(count ~ area2 + date3, data = eq2, family = poisson())
summary(eq_glm2)
summary(eq_glm)
plot(eq2$date3[eq2$area2=="Kantou"], eq2$count[eq2$area2=="Kantou"], type="l")
plot(eq2$date3[eq2$area2=="Kanto"], eq2$count[eq2$area2=="Kanto"], type="l")
plot(eq2$date3[eq2$area2=="Touhoku"], eq2$count[eq2$area2=="Touhoku"], type="l")

いつも奇麗なプログラムをアップしていますが、実際は変数名を間違えたりタイプミスがあったり、エラーを起こしながらプログラミングしています。そのうちエラー勉強会として題材にしても良いかも。。。？