このページは公開されるので個人情報などは記載しないこと。
課題が完成したらknitしPublishする。

課題名

rm(list = ls()) # 全オブジェクト削除
# ここにRコードを記入する。
#例）
plot(1:9)

# PCR移動平均

d <- read.csv('https://www.mhlw.go.jp/content/001060467.csv')
library(DT)

## Warning: パッケージ 'DT' はバージョン 4.3.2 の R の下で造られました

datatable(d,caption = 'PCR検査結果', colnames = c('日付', 'PCR検査実施件数'))

x <- as.POSIXct(d$日付, format = '%Y/%m/%d')
y <- d[,2]
yhat <- zoo::rollmean(y, 31, na.pad = T)
matplot(x,y, type = 'l', col = 3, 
        main = '新規陽性者数の推移（日別）', xlab = '年', ylab = '陽性者数')
matlines(x = x, y = yhat, col = 1)
grid()
legend('topleft', col = c(3, 1), lty = 1, legend = c('原系列', '31日移動平均'))

library(quantmod)

## Warning: パッケージ 'quantmod' はバージョン 4.3.2 の R の下で造られました

##  要求されたパッケージ xts をロード中です

## Warning: パッケージ 'xts' はバージョン 4.3.2 の R の下で造られました

##  要求されたパッケージ zoo をロード中です

## Warning: パッケージ 'zoo' はバージョン 4.3.2 の R の下で造られました

## 
##  次のパッケージを付け加えます: 'zoo'

##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     as.Date, as.Date.numeric

##  要求されたパッケージ TTR をロード中です

## Warning: パッケージ 'TTR' はバージョン 4.3.2 の R の下で造られました

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

z <- getSymbols('3092',src='yahooj', from = '2022-01-01',to   = '2022-12-31',auto.assign = F)
write.zoo(z, file = 'ZOZO.csv', sep = ',', quote = F)

x <- as.POSIXct(index(z), format = '%Y-%m-%d')
y <- z[, 'YJ3092.Close']
ylag1 <- lag(y, k = 1)
yhat5   <- rollmean(ylag1,  5, fill = NA, align = 'right')
yhat25  <- rollmean(ylag1, 25, fill = NA, align = 'right')
yhat25c <- rollmean(y, 25, fill = NA)
matplot(x, y, type = 'l', col = 3,
main = 'ZOZO（3092）', xlab = '日', ylab = '円')
matlines(x, yhat5,   col = 2)
matlines(x, yhat25,  col = 4)
matlines(x, yhat25c, col = 1)
grid()
legend('topleft', col = c(3, 2, 4, 1), lty = 1,
       legend = c('株価',
                  '5日移動平均',
                  '25日移動平均',
                  '25日中心化移動平均（参考）'))

# 予測精度

options(digits = 2)
set.seed(1)

n <- 10
x <- 1:n
y <- round(1:n + rnorm(n, mean = 1, sd = 2.5), 2)
e <- rnorm(n, mean = 0, sd = 1)
yhat <- yhat.out <- round(y + e, 2)
yhat.out[8] <- yhat.out[8] - 10 # 大外し

d <- data.frame(y, yhat, yhat.out)

library(DT)
datatable(d, colnames = c('観測値', '予測値', '予測値（大外し含む）'))

matplot(x = x, y = y, type = 'n', ylim = range(c(y, yhat, yhat.out)))
matlines(x = x, y = y, type = 'o', pch = 1, col = 1)
matlines(x = x, y = yhat, type = 'o', pch = 2, col = 2)
matlines(x = x, y = yhat.out, type = 'o', pch = 3, col = 3)

legend('topleft', pch = 1:3, col = 1:3, 
       legend = c('観測値', '予測値', '予測値（大外し含む）'))

get.accuracy <- function(yhat, y)
{
  data.frame(MBE  = mean(yhat - y),
             MAE  = mean(abs(yhat - y)),
             MAPE = mean(abs((yhat - y) / y)) * 100,
             RMSE = sqrt(mean((yhat - y)^2)))
}


(a <- get.accuracy(yhat, y))

##    MBE  MAE MAPE RMSE
## 1 0.25 0.83   46    1

(a.out <- get.accuracy(yhat.out, y))

##     MBE MAE MAPE RMSE
## 1 -0.75 1.6   53    3

a.out / a

##   MBE MAE MAPE RMSE
## 1  -3   2  1.2  2.9

a$MBE

## [1] 0.25

高め予測

機械学習

演習課題

2023-12-11

課題名