はじめに
先日、競馬予測関数を作成したのですが、予測と実際の結果を比較し、考察した結果、改善点が多くあるのではないかと考えられました。
—Rで競馬予測関数を作って実際に使用してみた【ウキゴリラの趣味】
そこで、改善点を修正した新しい競馬予測関数を作成し、検証した結果、有用である可能性が示唆されました！
—Rで作成した競馬予測関数がどのくらい有用なのか検証してみた【ウキゴリラの趣味】
本記事ではその競馬予測関数のプログラムコードについて公開したいと思います。
モデルの概要
本モデルでは重回帰分析を利用し、過去の半年間のレース結果から、走破タイムを予測する関数です。
応答変数に走破タイム
説明変数にレース距離、馬場状態、斤量、対象レースに出場する各馬、各馬とレース距離の相互作用（各馬に距離の得意不得意があるかを示す）を入れています。
このモデル式を使用して、対象レースの各馬の走破タイムを予測するモデルとなっています。
また、モデルの当てはまりの良さを示す赤池情報量基準（AIC）によるモデル選択を行うと、各馬とレース距離の相互作用は必要ないと判定されることが多いです。
そこで、説明変数にレース距離、馬場状態、斤量、対象レースに出場する各馬のみをいれたモデルでも予測を行います。
使い方
・使用するデータフレーム
以下の記事に公開した関数を用いて、各馬のデータを吸い出します。
—地方競馬のデータをRを使ったWebスクレイピングによって取得する方法【ウキゴリラの趣味】
“rbind”を用いて、そのレースで出場する全ての馬のデータを繋げます。
その際、デフォルトでは馬番(ID)は1となっているため、各馬に合わせて変換していきます。
自分は、以下のコードを用いて、一頭ずつデータを繋げています。
keibatyusyutu("対象の馬のURL")
dat<-subset(dat1_1,place=="対象の競馬場") #同じ競馬場のデータのみ抽出
dat$ID<-rep(1,length(dat$ID)) #馬番を変更
data<-rbind(data,dat) #全てのデータを繋げるまた、当日の斤量、馬番、レース距離、馬場状態を以下のようなデータフレームでまとめます。
この列名は変更するとエラーが生じるため、ご注意ください。
Todaydata<-data.frame(
Weight=c(56,51,53,56,56,53,56,56,56,56),
ID=c(1,2,3,4,5,6,7,8,9,10),
Dis=rep(1400,10),
Baba=rep("不",10)
)・関数の入力
上記でまとめたデータフレームを使い、以下のコードを入力します。
なお、各馬の過去のレース結果のデータをまとめたデータフレームは“data”、当日の情報をまとめたデータフレームを“Todaydata“としております。
Keiba.Time.Predict(data$time,data$dis,data$Weight,data$Baba,data$ID,Todaydata,"SAGA_040222_10R")
#Keiba.Time.Predict(走破タイム,レース距離,斤量,馬場状態,馬番,当日の情報をまとめたデータフレーム,レース名)の順に打ち込むすると、以下のような図が表示されると思います。
総合と書かれた図は、各馬の距離に応じた走破タイムの違い（距離と各馬の相互作用）を考慮していないモデルでの予測です。
これまでのレース結果から、馬場や斤量、レース距離の影響を統計的に排除したうえで、どの馬が一番早かったのかを比較することができます。
距離が書かれた図は、各馬の距離に応じた走破タイムの違い（距離と各馬の相互作用）を考慮したモデルでの予測です。
各馬は距離が長い方が強いのか、短い方が強いのか、検証する際にご活用ください。
さらに、“Result”と入力することによって、各馬の予測した走破タイムも出力することができます。
> Result
ID Time800 Time1300 Time1400 Time1600 Time1800 Time2000 TimeSogo
1 1 44.71224 84.55232 92.52034 108.4564 124.3924 140.3284 100.43641
2 2 56.63216 86.98866 93.05996 105.2026 117.3452 129.4878 100.76205
3 3 47.77942 85.62204 93.19057 108.3276 123.4647 138.6017 100.10018
4 4 50.81363 84.82998 91.63325 105.2398 118.8463 132.4529 99.63700
5 5 48.83251 84.95996 92.18545 106.6364 121.0874 135.5384 100.02207
6 6 50.87204 84.37204 91.07204 104.4720 117.8720 131.2720 98.43616
7 7 50.29115 85.11989 92.08564 106.0171 119.9486 133.8801 99.90362
8 8 55.38074 86.62998 92.87983 105.3795 117.8792 130.3789 101.45364
9 9 52.00384 86.00995 92.81117 106.4136 120.0161 133.6185 101.04875
10 10 50.43261 85.27494 92.24340 106.1803 120.1173 134.0542 100.27524プログラムコード
本モデルのプログラムコードは以下の通りです。
各々、自分好みにアレンジして、ご活用ください！
私の記事が皆様の力添えになれば幸いです。
Keiba.Time.Predict<-function(time,distance,Weight,Baba,IDs,todaydata,racename)
{
da<-data.frame(
Dis=sub("ダ", "", distance),
Weight=Weight,
Baba=Baba,
ID=IDs,
SecTime=
(as.numeric(as.POSIXlt(time, format="%M:%S", tz = "Japan")$min*60)
+as.numeric(as.POSIXlt(time, format="%M:%S", tz = "Japan")$sec)
+as.numeric(substr(time, 6, 6))/10)
)
predictdata<-data.frame(
Dis=rep(800,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
da<-na.omit(da)
Res<<-lm(SecTime~as.numeric(Dis)+as.numeric(Weight)+Baba+as.factor(ID)+as.numeric(Dis)*as.factor(ID),data=da)
predictdata<-data.frame(
Dis=rep(800,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
Time800<<-predict(Res
,predictdata
,interval="confidence"
, level=0.95
)
predictdata<-data.frame(
Dis=rep(1300,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
Time1300<<-predict(Res
,predictdata
,interval="confidence"
, level=0.95
)
predictdata<-data.frame(
Dis=rep(1400,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
Time1400<<-predict(Res
,predictdata
,interval="confidence"
, level=0.95
)
predictdata<-data.frame(
Dis=rep(1600,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
Time1600<<-predict(Res
,predictdata
,interval="confidence"
, level=0.95
)
predictdata<-data.frame(
Dis=rep(1800,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
Time1800<<-predict(Res
,predictdata
,interval="confidence"
, level=0.95
)
predictdata<-data.frame(
Dis=rep(2000,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
Time2000<<-predict(Res
,predictdata
,interval="confidence"
, level=0.95
)
Res2<<-lm(SecTime~as.numeric(Dis)+as.numeric(Weight)+Baba+as.factor(ID),data=da)
predictdata<-data.frame(
Dis=rep(1500,length(todaydata$ID)),
Weight=todaydata$Weight,
ID=todaydata$ID,
Baba=todaydata$Baba)
TimeSOGO<<-predict(Res2
,predictdata
,interval="confidence"
, level=0.95
)
Result<<-data.frame(
ID=todaydata$ID,
Time800=Time800[,1],
Time1300=Time1300[,1],
Time1400=Time1400[,1],
Time1600=Time1600[,1],
Time1800=Time1800[,1],
Time2000=Time2000[,1],
TimeSogo=TimeSOGO[,1])
plot(Result$ID,Result$Time800,xlab="馬番",
ylab="走破タイム",main=c(racename,"800m"),pch=19,
ylim=c((min(Time800[,2])-1),(max(Time800[,3])+1)))
arrows(todaydata$ID,Time800[,2],
todaydata$ID,Time800[,3],angle=90,code=3,length=0.1)
plot(Result$ID,Result$Time1300,xlab="馬番",
ylab="走破タイム",main=c(racename,"1300m"),pch=19,
ylim=c((min(Time1300[,2])-1),(max(Time1300[,3])+1)))
arrows(todaydata$ID,Time1300[,2],
todaydata$ID,Time1300[,3],angle=90,code=3,length=0.1)
plot(Result$ID,Result$Time1400,xlab="馬番",
ylab="走破タイム",main=c(racename,"1400m"),pch=19,
ylim=c((min(Time1400[,2])-1),(max(Time1400[,3])+1)))
arrows(todaydata$ID,Time1400[,2],
todaydata$ID,Time1400[,3],angle=90,code=3,length=0.1)
plot(Result$ID,Result$Time1600,xlab="馬番",
ylab="走破タイム",main=c(racename,"1600m"),pch=19,
ylim=c((min(Time1600[,2])-1),(max(Time1600[,3])+1)))
arrows(todaydata$ID,Time1600[,2],
todaydata$ID,Time1600[,3],angle=90,code=3,length=0.1)
plot(Result$ID,Result$Time1800,xlab="馬番",
ylab="走破タイム",main=c(racename,"1800m"),pch=19,
ylim=c((min(Time1800[,2])-1),(max(Time1800[,3])+1)))
arrows(todaydata$ID,Time1800[,2],
todaydata$ID,Time1800[,3],angle=90,code=3,length=0.1)
plot(Result$ID,Result$TimeSogo,xlab="馬番",
ylab="走破タイム",main=c(racename,"総合"),pch=19,
ylim=c((min(TimeSOGO[,2])-1),(max(TimeSOGO[,3])+1)))
arrows(todaydata$ID,TimeSOGO[,2],
todaydata$ID,TimeSOGO[,3],angle=90,code=3,length=0.1)
da<<-da
}
コメント