modelisation_predictive/_04_aggregation.R at master · raphbrd/modelisation_predictive · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#' Fit agregation of previous models
#' Agregation of a boosting model and three gam models

#1st expert
equation <- Net_demand~  Time + toy + Temp + Net_demand.1 + Net_demand.7 + Temp_s99 + as.factor(WeekDays3) +
            Wind_power.1 + Solar_power.1 + Wind_power.7 + Solar_power.7 + Wind_weighted + Nebulosity_weighted +
            confinement1 + confinement2 + confinement3

Ntree <- 1000
gbm1 <- gbm(equation,distribution = "gaussian" , data = Data0[sel_a,],n.trees = Ntree, interaction.depth = 10,
            n.minobsinnode = 5, shrinkage = 0.02, bag.fraction = 0.5, train.fraction = 1,
            keep.data = FALSE, n.cores = 8)

best.iter=gbm.perf(gbm1,method="OOB", plot.it = TRUE,oobag.curve = TRUE)
best.iter
NTreeOpt <- which.min(-cumsum(gbm1$oobag.improve))
gbm1.forecast <- predict(gbm1,n.trees=NTreeOpt,single.tree=FALSE,newdata=Data0[sel_b,])

#2nd expert
equation <- Net_demand ~ s(as.numeric(Date),k=3, bs='cr') + s(toy,k=30, bs='cc') +
  s(Temp,k=10, bs='cr') + s(Load.1, bs='cr')+ s(Load.7, bs='cr')  +
  s(Temp_s99,k=10, bs='cr') + WeekDays + s(Wind_weighted) +
  te(as.numeric(Date), Nebulosity, k=c(4,10)) +
  s(Wind_power.1, k=10, bs='cr') + s(Solar_power.1, k=10,  bs='cr')


md.gam.1 <- gam(equation,  data=Data0[sel_a,])
md.gam.1.forecast <- predict(md.gam.1,  newdata = Data0[sel_b, ])

#3rd expert
equation = Net_demand ~ s(as.numeric(Date), k = 3, bs = 'cr') +
  s(toy, k = 30, bs = 'cc') + s(Temp, k = 10, bs = 'cr') +
  s(Net_demand.1, bs = 'cr',  by = WeekDays) +
  s(Net_demand.7, bs = 'cr') +
  s(Load.1, bs = 'cr',  by = WeekDays) +
  s(Load.7, bs = 'cr') +
  s(Temp_s99, k = 10, bs = 'cr') +
  as.factor(WeekDays) +
  as.factor(BH) +
  s(Wind_weighted) +
  te(as.numeric(Date), Nebulosity_weighted, k = c(4, 10)) +
  s(Wind_power.1, k = 10, bs = 'cr') +
  s(Solar_power.1, k = 10,  bs = 'cr')

md.gam.2 <- gam(equation,  data=Data0[sel_a,])
md.gam.2.forecast <- predict(md.gam.2,  newdata = Data0[sel_b, ])

#4th expert
equation = Net_demand ~ s(as.numeric(Date), k = 3, bs = 'cr') +
  s(toy, k = 30, bs = 'cc') + s(Temp, k = 10, bs = 'cr') +
  s(Net_demand.1, bs = 'cr',  by = WeekDays) +
  s(Net_demand.7, bs = 'cr') +
  s(Load.1, bs = 'cr',  by = WeekDays) +
  s(Load.7, bs = 'cr') +
  s(Temp_s99, k = 10, bs = 'cr') +
  WeekDays +
  s(Wind_weighted) +
  te(as.numeric(Date), Nebulosity_weighted, k = c(4, 10))

md.gam.3 <- gam(equation,  data=Data0[sel_a,])
md.gam.3.forecast <- predict(md.gam.3,  newdata = Data0[sel_b, ])

# Agregation

experts = cbind(gbm1.forecast,md.gam.1.forecast,md.gam.2.forecast,md.gam.3.forecast)
colnames(experts) = c("gbm","gam1","gam2","gam3")
or = oracle(Data0$Net_demand[sel_b],experts,model="convex",loss.type="square")
pinball_loss(Data0$Net_demand[sel_b],or$prediction,quant = 0.95)

# Return a pinball loss of 107, the same score as our best gam model on the validation test, all the weigth is on this model (gam2)

# Agregation removing the best gam and the boosting model and whithout quantile prediction

equation1 <- Net_demand ~ s(as.numeric(Date),k=3, bs='cr') + s(toy,k=30, bs='cc') +
  s(Temp,k=10, bs='cr') + s(Load.1, bs='cr')+ s(Load.7, bs='cr')  +
  s(Temp_s99,k=10, bs='cr') + WeekDays + s(Wind_weighted) +
  te(as.numeric(Date), Nebulosity, k=c(4,10)) +
  s(Wind_power.1, k=10, bs='cr') + s(Solar_power.1, k=10,  bs='cr')


md.gam.1 <- gam(equation1,  data=Data0[sel_a,])
md.gam.1.forecast <- predict(md.gam.1,  newdata = Data0[sel_b, ])


equation3 = Net_demand ~ s(as.numeric(Date), k = 3, bs = 'cr') +
  s(toy, k = 30, bs = 'cc') + s(Temp, k = 10, bs = 'cr') +
  s(Net_demand.1, bs = 'cr',  by = WeekDays) +
  s(Net_demand.7, bs = 'cr') +
  s(Load.1, bs = 'cr',  by = WeekDays) +
  s(Load.7, bs = 'cr') +
  s(Temp_s99, k = 10, bs = 'cr') +
  WeekDays +
  s(Wind_weighted) +
  te(as.numeric(Date), Nebulosity_weighted, k = c(4, 10))

md.gam.3 <- gam(equation3,  data=Data0[sel_a,])
md.gam.3.forecast <- predict(md.gam.3,  newdata = Data0[sel_b, ])


experts = cbind(md.gam.1.forecast,md.gam.3.forecast)
colnames(experts) = c("gam1","gam3")
or = oracle(Data0$Net_demand[sel_b],experts,model="convex",loss.type="square")
pinball_loss(Data0$Net_demand[sel_b],or$prediction,quant = 0.95)

# Return a pinball loss 115 on the validation test, 0.1 weight on gam1 and 0.9 on gam3
# Using this agregation model to predict the private set data1

md.gam.1.data1 = gam(equation1,  data=Data0)
md.gam.1.data1.forecast = predict(md.gam.1.data1, newdata=Data1)

md.gam.3.data1 = gam(equation3,  data=Data0)
md.gam.3.data1.forecast = predict(md.gam.3.data1, newdata=Data1)

experts = cbind(md.gam.1.data1.forecast,md.gam.3.data1.forecast)

agreg_gam1_gam3 = predict(or,newexpert = experts)

# This model generalize well and gives us good score on both public and private Kaggle's datasets