Termékek forgalma viharos időben a Walmartnál

Pethő Gergely

1. Bevezetés

Dolgozatomban a Kaggle 2015. májusi versenyfeladatát próbáltam megoldani. A feladatot a Walmart írta ki adatelemzők toborzása céljából.

A feladat lényege, hogy meg kell próbálni előrejelezni bizonyos termékek forgalmát az időjárás függvényében. A feladat kiírása szerint olyan termékek tartoznak az időjárástól potenciálisan függő forgalmú áruk közé, mint az esernyő (értelemszerűen), illetve a tej és a kenyér. Várható, hogy viharok idején megnő például az esernyők forgalma, az áruházak készletezői számára viszont problémát jelent, hogy megbecsüljék, pontosan mekkora mennyiségre van szükség ilyenkor az adott áruból, hogy se ne fogyjanak ki idő előtt a készletből, se ne rendeljenek túl sokat.

A feladat a következő oldalon érhető el: Walmart Recruiting II: Sales in Stormy Weather

A feladat megoldásához felhasználható adatok négy CSV-fájlban lettek közzéadva a Kaggle oldalán. Ezek a fájlok a következők: * weather.csv: egy aránylag sűrű mátrix, amely 20 meteorológiai állomás által mért időjárási adatokat tartalmazza a vizsgált időszak minden napján (2012. januártól 2014. októberig, ez durván 1000 nap): a mérés dátumát, hőmérsékleti adatokat, a napkelte és a napnyugta idejét, az adott napi csapadék típusát, a hó, illetve az összes csapadék mennyiségét, légnyomást, valamint a szél irányát és sebességét. Mérete kb. 2 MB. * train.csv: egy viszonylag ritka mátrix egy dimenzióba kilapított változata, amelynek soraiban az áll, hogy egy adott napon a vizsgált 45 Walmart-áruház egyikében a vizsgált 111 árucikk egyikéből hány darabot adtak el. Az adathalmaz leírása megjegyzi, hogy különböző áruházakban vagy térségekben ugyanaz az árucikk más-más kóddal szerepelhet; és ugyan ezt külön nem teszi hozzá, de sejthető, hogy ez fordítva is igaz lehet, tehát előfordulhat, hogy különböző áruházakban ugyanaz a kód különböző árucikket jelenthet. A fájl mérete kb. 80 MB. * key.csv: egy két oszlopból álló kapcsolótábla, amely megadja, hogy a train.csv táblázatban szereplő áruházak egyenként melyik meteorológiai állomáshoz tartoznak a weather.csv táblázatból. * test.csv: formátuma hasonló a train.csv-hez, tehát minden sorában egy dátum, egy áruház és egy árucikk sorszáma szerepel, viszont hiányzik az adott árucikk eladott mennyisége. A feladat ennek az adatnak a megbecslése.

Pontosabban a feladatot úgy határozták meg, hogy a vizsgált időszak nagyjából első felében (2013 áprilisig) a “fontos időjárási események” 3 napos környezetébe eső eladási adatok a tanító adathalmazban szerepelnek, míg az időszak második felében az ilyen napok eladási adatai hiányoznak, és ezek alkotják a teszthalmazt, jobban mondva ezek az előrejelzendő adatok. A “fontos időjárási esemény” olyan napot jelent, amelyen legalább 1 hüvelyk hó vagy eső esett.

Fontos megjegyezni, hogy mind az üzletek, mind a meteorológiai állomások, mind a termékek azonosítása pusztán kódszámmal történik a feladatban, amelyhez a Walmart nem adott kulcsot. Ez azt jelenti, hogy a létesítmények tényleges földrajzi helyzetéről nem tudunk semmit, csupán a key.csv fájl mondja meg nekünk, hogy melyik üzletek esnek egy-egy állomás körzetébe. A termékeknek pedig pusztán a kódszámát és az eladásaik számát ismerjük, azt viszont nem tudjuk, hogy a kódok milyen jellegű termékeket fednek, például esernyőt vagy éppenséggel naptejet.

A leglényegesebb fájl, tehát a train.csv olyan nagy, hogy a RapidMiner 1 GB-os korlátozott memóriájába nem fér bele. Részben ezért úgy döntöttem, hogy a dolgozathoz az R-t használom, és a dolgozatot az RStudioban rendelkezésre álló R Markdownban írom meg.

2. Az adatok előkészítése

2.1. train.csv

Mindenekelőtt feltűnő, hogy a train.csv fájlban az adatok igen pazarló módon vannak megadva. A táblázat a (dátum, áruház, árucikk) hármasokhoz akkor is tartalmaz egy-egy sort, ha a sorban történetesen 0 érték szerepel az eladási oszlopban. Ha betöltjük a táblázatot, és megszámoljuk a 0-nál nagyobb darabszámú értékesítéseket, kiderül, hogy a fájl több mint 4,5 millió sorából valójában csak kb. 120 ezerben van nem redundáns adat:

t <- read.csv('f:/downloads/train.csv',colClasses=c("Date","factor","factor","integer"))
str(t)
## 'data.frame':    4617600 obs. of  4 variables:
##  $ date     : Date, format: "2012-01-01" "2012-01-01" ...
##  $ store_nbr: Factor w/ 45 levels "1","10","11",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ item_nbr : Factor w/ 111 levels "1","10","100",..: 1 24 35 46 57 68 79 90 101 2 ...
##  $ units    : int  0 0 0 0 0 0 0 0 29 0 ...
head(t)
##         date store_nbr item_nbr units
## 1 2012-01-01         1        1     0
## 2 2012-01-01         1        2     0
## 3 2012-01-01         1        3     0
## 4 2012-01-01         1        4     0
## 5 2012-01-01         1        5     0
## 6 2012-01-01         1        6     0
nrow(t[t[4] > 0,])
## [1] 118696

Ha az adatokat a jelenlegi (ún. hosszú) adatszerkezetből átrendezzük egy olyanba, ahol csupán a bolt-dátum párokra esik egy-egy sor, és a termékeket oszlopokba tesszük, akkor a 4,5 millió sor helyett kb. 41 ezret kapunk. Bár az oszlopok száma 4-ről 113-ra nő, ez még így is jelentős megtakarítást jelent, a memóriaigény várhatóan negyedére esik vissza, ugyanakkor bizonyos függvények ezzel a (széles) táblázattal nem tudnak dolgozni. Ennek ellenére érdemes elvégezni ezt az átalakítást, mert áttekinthetőbbé válik az adathalmaz. Ugyanakkor így CSV-be mentve már csak tizedakkora az adatfájl, mint eredetileg volt, így már betölthető RapidMinerbe is.

library(reshape2)
szeles <- dcast(t,date + store_nbr ~ item_nbr)
## Using units as value column: use value.var to override.
names(szeles) <- c("date","store_nbr",names(szeles[3:113]))

Tekintettel arra, hogy a boltok egészen különböző területeken fekszenek (és így eltérő az időjárásuk), továbbá az árukódszámaik sem feltétlenül egyeznek, az egyik boltra vonatkozó eladási adatok várhatóan nem hordoznak információt arra vonatkozóan, hogy egy másik boltban milyenek lesznek az eladások. Emiatt a legésszerűbb minden egyes bolt adatsorát külön táblázatként kezelni a későbbiekben. Ez a következő módon oldható meg, példaként az 1. számú üzletet véve:

s <- szeles[szeles$store_nbr=="1",]
datetemp <- s[1]
s <- s[3:113]
s <- as.matrix(s)

Megvizsgálhatjuk például, hogy hányféle árucikkből adtak el egyáltalán az adott üzletben, és egyenként összesen mennyit.

s <- s[,colSums(s) != 0]
ncol(s)
## [1] 8
colSums(s)
##    28    40    47    51    89     9    93    99 
##  4893   254  2409   925   157 27396  1103   492

Áttekinthetjük az egyes termékek idősorait is. Azt láthatjuk például, hogy a 28. termék időben viszonylag stabil, nincsenek feltűnő szezonális eltérések.

s <- cbind(datetemp,s)
plot(s$date,s$"28", type="l",main="A 28. termék forgalma",xlab="dátum",ylab="eladott mennyiség")

A 40. termék eladásai láthatóan őszre és télre koncentrálódnak, egyértelműen szezonálisak. Várhatóan az ilyenek fognak leginkább függni az időjárási eseményektől.

plot(s$date,s$"40", type="l",main="A 40. termék forgalma",xlab="dátum",ylab="eladott mennyiség")

Speciális problémát jelez a 93. termék. Ebből egész évben adtak el 2012-ben, nyáron valamivel többet, mint télen, viszont 2013 elejétől egy darabot sem, ami arra utal, hogy az adott üzletben egyáltalán nem forgalmazták már ezt az árucikket, talán meg is szűnt. Az előrejelzéseinknél ezt figyelembe kell majd venni.

plot(s$date,s$"93", type="l",main="A 93. termék forgalma",xlab="dátum",ylab="eladott mennyiség")

Épp fordított problémát jelez a 99. termék. Ebből 2012 közepén eladtak néhány darabot, de egyébként a vizsgált időszak első felében semmit, utána 2013 őszétől viszont különösebb ingadozás nélkül folyamatosan forgalmazták. Ez azért okoz nehézséget, mert a 2013 április utáni viharos napokra vonatkozó előrejelzéseinket a feladat értelmében részben a 2013 áprilisig terjedő viharos napok eladásai alapján kellene megfogalmaznunk, az ehhez hasonló esetekben azonban az utóbbiak vélhetően 0-ra jönnek ki, míg a tényleges érték vélhetően inkább a forgalmazási időszak (2013 őszétől 2014 végéig) egy átlagos napjának a forgalmához fog közel állni.

plot(s$date,s$"99", type="l",main="A 99. termék forgalma",xlab="dátum",ylab="eladott mennyiség")

2.2. key.csv

A következő lépés a kapcsolótábla betöltése, ez speciális feldolgozást nem igényel.

key <- read.csv('f:/downloads/key.csv',colClasses=c("factor","factor"))
str(key)
## 'data.frame':    45 obs. of  2 variables:
##  $ store_nbr  : Factor w/ 45 levels "1","10","11",..: 1 12 23 34 41 42 43 44 45 2 ...
##  $ station_nbr: Factor w/ 20 levels "1","10","11",..: 1 6 18 20 4 6 17 15 9 4 ...

Áttekinthetjük, hogy melyik meteorológiai állomáshoz hány Walmart-üzlet tartozik az adathalmazban:

a <- aggregate(key$store_nbr,list(key$station_nbr),length)
names(a)[1] <- "station_nbr"
names(a)[2] <- "number_of_shops"
t(a)
##                 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## station_nbr     "1"  "10" "11" "12" "13" "14" "15" "16" "17" "18"  "19" 
## number_of_shops "1"  "3"  "2"  "4"  "5"  "4"  "1"  "2"  "6"  "1"   "1"  
##                 [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## station_nbr     "2"   "20"  "3"   "4"   "5"   "6"   "7"   "8"   "9"  
## number_of_shops "1"   "1"   "3"   "1"   "1"   "2"   "3"   "1"   "2"
mean(a$number_of_shops)
## [1] 2.25

Láthatjuk, hogy az üzletek száma állomásonként 1 és 6 közé esik (a 17. állomás esetében 6), és átlagosan nagyjából 2 üzlet tartozik egy állomáshoz.

2.3. weather.csv

Az időjárási adatok betöltése ezzel szemben komolyabb előfeldolgozást követel meg. A nyers adathalmazban számos olyan mező (oszlop) található, amelyek a feladat szempontjából teljesen irrelevánsnak tűnnek. Legalább a következő mezők ezért kizárhatóak: a hőmérséklet normálistól való eltérése, a harmatpont, a “nedves hőmérséklet” (wet bulb), a melegedés, hűlés, napkelte, napnyugta, a légnyomás és a szélirány. A hőmérsékleti értékeket (napi minimum, maximum, átlag), az átlagos szélsebességet és természetesen a csapadékadatokat (fajtája és mennyisége) ellenben célszerű megtartani, mert ezek valószínűleg befolyásolhatják a vásárlói magatartást az időjárási eseményekkel összefüggésben.

A számadatokat sztringként olvasom be a fájlból. Ennek két oka van: egyrészt idézőjelekkel vannak a fájlban körülvéve, ezért az R nem hajlandó ezeket eleve számként értelmezni (tehát hibaüzenetet ad, ha arra utasítom, hogy numerikus értékként olvassa be ezeket), másrészt a csapadékmennyiség oszlopban az elhanyagolható mennyiséget a “T” érték jelöli. A hőmérsékleti (3-5. oszlop) és szélerősségre vonatkozó (9. oszlop) értékeket a sztringként való beolvasás után azonnal átalakítom számokká. A hó, illetve összes csapadék mennyiségére vonatkozó adatokból az átalakítás előtt számmá írom át a “T” értékeket, ezeket 0-nak tekintem.

weather <- read.csv('f:/downloads/weather.csv', colClasses=c("factor", "Date", "character", "character", "character", "NULL","NULL","NULL","NULL","NULL","NULL","NULL", "character", "character", "character", "NULL","NULL","NULL","NULL","character"), na.strings=c("M"),strip.white=TRUE)
weather[3] <- as.integer(unlist(weather[3]))
weather[4] <- as.integer(unlist(weather[4]))
weather[5] <- as.integer(unlist(weather[5]))
weather[9] <- as.numeric(unlist(weather[9]))
weather[grep("T",unlist(weather[7])),7] <- 0
weather[grep("T",unlist(weather[8])),8] <- 0
weather[7] <- as.numeric(unlist(weather[7]))
weather[8] <- as.numeric(unlist(weather[8]))
str(weather)
## 'data.frame':    20517 obs. of  9 variables:
##  $ station_nbr: Factor w/ 20 levels "1","10","11",..: 1 12 14 15 17 18 19 20 2 3 ...
##  $ date       : Date, format: "2012-01-01" "2012-01-01" ...
##  $ tmax       : int  52 48 55 63 63 50 66 34 73 72 ...
##  $ tmin       : int  31 33 34 47 34 33 45 19 53 48 ...
##  $ tavg       : int  42 41 45 55 49 42 NA 27 63 60 ...
##  $ codesum    : chr  "RA FZFG BR" "RA" " " " " ...
##  $ snowfall   : num  NA 0 0 0 0 0 NA NA NA 0 ...
##  $ preciptotal: num  0.05 0.07 0 0 0 0 0 0 0 0 ...
##  $ avgspeed   : num  4.6 11.3 10 8.2 13.8 10.2 10.9 22.5 5.5 4.8 ...

Azon napok száma, amikor a feladatban meghatározott értelemben nagy mennyiségű eső, illetve hó esett, a következő módon számítható ki. A 6. oszlopban az “RA” kód esőt, az “SN” havat jelent.

nrow(weather[grepl("RA",unlist(weather$codesum)) & weather$preciptotal >= 1 & !is.na(weather$preciptotal),])
## [1] 516
nrow(weather[grepl("SN",unlist(weather$codesum)) & weather$snowfall >= 1 & !is.na(weather$snowfall),])
## [1] 140

Az adathalmazt alaposabban megvizsgálva feltűnik, hogy míg minden más meteorológiai állomáshoz 1035 napnyi mérés tartozik, addig az 5. állomáshoz csupán 852 napnyi.

t(aggregate(weather,list(weather$station_nbr),length)[1:2])
##             [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  
## Group.1     "1"    "10"   "11"   "12"   "13"   "14"   "15"   "16"   "17"  
## station_nbr "1035" "1035" "1035" "1035" "1035" "1035" "1035" "1035" "1035"
##             [,10]  [,11]  [,12]  [,13]  [,14]  [,15]  [,16]  [,17]  [,18] 
## Group.1     "18"   "19"   "2"    "20"   "3"    "4"    "5"    "6"    "7"   
## station_nbr "1035" "1035" "1035" "1035" "1035" "1035" " 852" "1035" "1035"
##             [,19]  [,20] 
## Group.1     "8"    "9"   
## station_nbr "1035" "1035"

Ezt történetesen leolvasható a Kaggle-oldalon a feladat adathalmazát dokumentáló lapon is. Az ott közölt ábrán látszik, hogy az 5. állomáshoz hiányzik az időjárási adatok egy része, viszont cserébe ehhez az állomáshoz nem is tartoznak olyan napok, amelyeken előre kellene jelezni az időjárást (ezeket pirossal jelöli az ábra).

library(png)
pngfile <- "f:/downloads/weather_events.png"
pic <- readPNG(pngfile)
grid::grid.raster(pic)

Mindez gyakorlatilag azt jelenti, hogy az 5. állomáshoz tartozó üzletekhez nem szükséges előrejelzéseket készítenünk. Mivel az üzleteket egyébként is egyenként fogom vizsgálni a fent kifejtettek miatt, ezért az 5. állomáshoz tartozó üzletek adatai teljesen törölhetőek az adathalmazból.

key[key$station_nbr==5,]
##    store_nbr station_nbr
## 35        35           5

Láthatjuk, hogy egyetlen ilyen üzlet van. Kitörlöm az adathalmazból mind az 5. állomás, mind a 35. üzlet sorait.

weather <- weather[weather$station_nbr != 5,]
szeles <- szeles[szeles$store_nbr != 35,]

Érdemes egy pillantást vetni a meteorológiai állomások méréseinek összesített értékeire is.

allomasok <- aggregate(weather[c(5,9)],list(weather$station_nbr),mean,na.rm=TRUE)
temp <- aggregate(weather[7:8],list(weather$station_nbr),sum,na.rm=TRUE)
allomasok <- cbind(allomasok,temp)[c(1,2,3,5,6)]
names(allomasok)[1] <- "station_nbr"
allomasok
##    station_nbr     tavg  avgspeed snowfall preciptotal
## 1            1 52.30321  6.572190      0.0       98.64
## 2           10 73.11295  8.401741      0.0      153.27
## 3           11 69.70126  6.327422      0.0      146.39
## 4           12 71.36566  6.376573      0.0      170.69
## 5           13 51.76983  7.905319      0.0       46.32
## 6           14 62.97959 11.173959     19.7      107.21
## 7           15 48.16796  5.327422    186.1      116.88
## 8           16 51.09785  7.354528      0.0      114.39
## 9           17 62.58835  6.265474      0.0       86.21
## 10          18 69.51115  9.434366      5.1       83.48
## 11          19 49.31256  9.310358    136.7       98.55
## 12           2 50.01271  9.721435    217.4      133.73
## 13          20 68.87938  7.987803      0.0       85.05
## 14           3 62.66148  9.544240     17.8       86.67
## 15           4 71.37171  5.141433      0.2      102.46
## 16           6 69.31293  7.928848      0.0       92.53
## 17           7 63.39638 11.814438      0.0       98.36
## 18           8 73.26808  5.444660      0.0       69.76
## 19           9 47.81890  8.267087      0.0       83.74

A középhőmérséklet a vizsgált kb. 3 éves időszakban 47 és 73 °F (azaz kb. 8 és 22 °C) között mozog az egyes állomásokon, míg a 3 év összes csapadéka (az adatok alapján a hómennyiséget leszámítva) 46 és 170 hüvelyk, azaz éves szinten 410 mm és 1510 mm között mozog. Összehasonlításként Magyarország középhőmérséklete 9 és 11 °C, csapadékmennyisége pedig tájanként 500 és 750 mm között mozog. A 19 állomásból csak 7 mért hóesést, viszont a 2. állomáson éves szinten kb. 2 méter hó hullott, ami elég soknak hangzik. Az adatokból arra következtethetünk, hogy a meteorológiai állomások (és így az üzletek) olyan erőteljesen különböző területeken fekszenek, amelyek között nem csupán időjárási, hanem éghajlati jellegű különbségek is vannak. Ez ismét amellett szól, hogy az áruházak eladási adatait egymástól elkülönítve vizsgáljuk, és ne próbáljunk rájuk egy közös modellt építeni.

2.4. test.csv

Ennek a fájlnak a szerkezete megegyezik a tanító adathalmazéval. Attól mindössze abban különbözik, hogy az árucikkekhez nincs megadva eladott darabszám (ugyanis ezt a változót kell előrejeleznünk). Az adatok betöltése ezért nagyon hasonlóan történik.

test <- read.csv('f:/downloads/test.csv',colClasses=c("Date","factor","factor"))
str(test)
## 'data.frame':    526917 obs. of  3 variables:
##  $ date     : Date, format: "2013-04-01" "2013-04-01" ...
##  $ store_nbr: Factor w/ 44 levels "1","10","11",..: 12 12 12 12 12 12 12 12 12 12 ...
##  $ item_nbr : Factor w/ 111 levels "1","10","100",..: 1 24 35 46 57 68 79 90 101 2 ...
head(test)
##         date store_nbr item_nbr
## 1 2013-04-01         2        1
## 2 2013-04-01         2        2
## 3 2013-04-01         2        3
## 4 2013-04-01         2        4
## 5 2013-04-01         2        5
## 6 2013-04-01         2        6

Az előrejelzett adatokat ugyan egy külön vektorban is tárolhatnánk, de kényelmesebb erre a célra egy egyelőre kinullázott oszlopot létrehozni ebben a táblázatban.

test[,"units"] <- 0

2.5. Boltok és meteorológiai állomások adatainak összefésülése

Ahhoz, hogy előrejelzéseket tudjunk megfogalmazni, célszerű egy közös táblázatba beletenni az eladási és az időjárási adatokat. Ehhez egy SQL jellegű összekapcsolás kínálkozik.

szeles <- merge(szeles,key)
szeles <- merge(szeles,weather,by=c("station_nbr","date"))

Ez az összekapcsolás bizonyos mértékű redundanciát okoz a táblázatban (láttuk, hogy átlagosan minden meteorológiai állomáshoz két-két üzlet tartozik, így az időjárási adatokat átlagosan kétszer megismételjük a táblázatban), de az egész meteorológiai táblázat mindössze néhány MB-ot foglal, így ezzel együtt lehet élni, másfelől viszont jelentősen könnyebb így dolgozni az adatokkal.

Megjegyzendő, hogy az előrejelzendő adatok táblázatával ezzel szemben nem érdemes ezen a ponton összekapcsolni az időjárási adatokat, az ugyanis rengeteg redundanciával járna, és a tesztadatsor elemei túlnyomó többségének az előrejelzéséhez nem lesz szükségünk időjárási adatokra. Itt tehát a hasonló összekapcsolásokat csak a konkrét előrejelzések készítésekor fogom végrehajtani, és csak akkor, ha szükség van rá.

Miután összekapcsoltuk az eladási és időjárási adatokat, könnyedén készíthetünk az idősorokról diagramokat, amelyek együttesen ábrázolják ezeket. Például lássuk az 1. üzletben a 28. termék eladási számait a napi átlagos hőmérséklet görbéjével (piros vonal), illetve a napi csapadék görbéjével (zöld vonal). Az ábrán kék karikák jelzik a különlegesként definiált időjárási jelenségeket. A diagramon csak 2013. március végéig jelenítem meg az adatokat, mert eladási adataink a különleges időjárási események által érintett napokról 2013. április után (a feladat értelmében) nincsenek. Az ábra megerősíti azt az előzetes gyanúnkat, hogy a 28. termék eladásai nem látszanak korrelálni az időjárással.

s <- szeles[szeles$store_nbr=="1",c(2,4:117,119:121)]
datatemp <- s[c(1,113:118)]
s <- s[2:112]
s <- as.matrix(s)
s <- s[,colSums(s) != 0]
s <- cbind(datatemp,s)
s <- s[s$date < "2013-04-01",]

plot(s$date,s$"28", type="l",main="A 28. termék és az időjárás",xlab="dátum",ylab="eladott mennyiség")
lines(s$date,s$tavg/5,col="red")
lines(s$date,s$preciptotal * 10,col="green")
for (i in 1:nrow(s)) {
    if (!is.na(s[i,"preciptotal"]) & s[i,"preciptotal"] > 1) {
        points(s[i,"date"],s[i,"preciptotal"] * 10, col="blue")
    }
}
legend("topleft",legend=c("forgalom","átl. hőmérséklet","csapadék"),col=c("black","red","green"),lty=c(1,1,1),border="black")

Ezzel szemben az 51. árucikk bizonyos összefüggést mutat a hőmérséklettel (szinte csak hideg időszakokban veszik), ellenben semmi összefüggés nem látható a csapadék mennyiségével.

plot(s$date,s$"51", type="l",main="Az 51. termék és az időjárás",xlab="dátum",ylab="eladott mennyiség")
lines(s$date,s$tavg/5,col="red")
lines(s$date,s$preciptotal * 10,col="green")
for (i in 1:nrow(s)) {
    if (!is.na(s[i,"preciptotal"]) & s[i,"preciptotal"] > 1) {
        points(s[i,"date"],s[i,"preciptotal"] * 10, col="blue")
    }
}
legend("topleft",legend=c("forgalom","átl. hőmérséklet","csapadék"),col=c("black","red","green"),lty=c(1,1,1),border="black")

3. Előrejelzések

A fent kifejtettek alapján sem hatékony, sem az eredmények minősége szempontjából célravezető nem lenne, ha elkezdenénk ész nélkül regressziós modelleket generálni és alkalmazni az adathalmazra. A tesztadatsor nagy részének a kitöltéséhez modellekre egyáltalán nem lesz szükség.

3.1. A boltokban nem forgalmazott áruk kezelése

Fent láttuk, hogy az egyes boltok csak néhány árucikket forgalmaznak, az 1. üzlet történetesen a 111 árucikkből csupán 8-at. Az egyes boltokban nem forgalmazott áruk eladásai nyilván a viharos napokon is 0 szinten lesznek, ezért a teszttáblázat e sorait érdemes az eredetileg beállított 0 értéken hagyni, és csak a fennmaradókhoz nyúlni.

Ennek érdekében a fent látott módon bejárom a boltokat, egyenként megállapítom, hogy mely árucikkeket forgalmazzák egyáltalán, majd kitöltöm a tesztadatok megfelelő sorait NA értékekkel. Ennek az eredménye az lesz, hogy “units” alatt a nem forgalmazott termékekre az eredetileg beállított 0 értékek fognak szerepelni, míg a forgalmazott termékekre definiálatlan érték jelzi, hogy ott még van teendőnk.

termekek_szama_boltonkent <- rep(0,35)

for (b in c(1:45)) {
    s <- szeles[szeles$store_nbr==as.character(b),c(4:114)]
    s <- s[,colSums(s) != 0]
    termekek_szama_boltonkent[b] <- ncol(s)
    
    for (i in 1:111) {
        if (as.character(i) %in% colnames(s)) {
            test[test$store_nbr == as.character(b)
                 & test$item_nbr == as.character(i),
                 "units"] <- NA
        }   
    }
}

Megjegyzendő, hogy természetesen fordítva is el lehetne járni, tehát NA-ra inicializálni az előrejelzés-vektort, és 0-ra átírni, ha nem forgalmazzák a terméket, de ez jóval több munka lenne, tekintve, hogy sokkal több a nem forgalmazott, mint a forgalmazott termék boltonként. A következő vektor mutatja, hogy melyik üzlet hány árut forgalmaz.

names(termekek_szama_boltonkent) <- as.character(1:45)
termekek_szama_boltonkent
##  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 
##  8  5  6  7  6  4  5  6  6  5  6  5  5  8  8  8  4  4  8  6  4  6  4  5  7 
## 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 
##  3  5  4  4  6  6  7  5  6  0  4  7  6  4  7  6  3  6  6  6

Vegyük észre, hogy a 35. boltnál 0 érték szerepel, ugyanis a hozzá tartozó adatokat töröltük a tanító táblázatból. Látjuk viszont, hogy 8-nál több árucikket egyik üzlet sem forgalmaz, így a tesztadatsor nagyon nagy része kinullázva maradt a fenti körben.

nrow(test)
## [1] 526917
nrow(test[is.na(test$units),])
## [1] 26168

Azaz a tesztadatsor 527000 sorából félmilliót ebben a lépésben már ki is töltöttünk. A továbbiakban ezzel a kisebb adathalmazzal dolgozom.

test.munka <- test[is.na(test$units),]
test <- test[!is.na(test$units),]

3.2. Sehol sem forgalmazott áruk kezelése

Ugyan ez viszonylag váratlan, de vannak olyan áruk a 111 között, amelyekből egyik üzletben sem adtak el a dokumentált időtartamban egy darabot sem. Az alábbi módon kiszámolom, hogy az egyes termékekből összesen hány darabot adtak el az összes üzletben (kivéve a 35.-ben) összesen.

c <- colSums(szeles[4:114])
summary(c)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0     476     766   40710    2242 1005000
c[c<100]
## 102  24  63  66 
##  31   0   0   0

Láthatóan három ilyen termék van, továbbá a 102. számú termékből összesen 31-et adtak el 3 év alatt (ebből is csupán kettőt 2013. április előtt), tehát ennek forgalma szintén valószínűleg 0 lesz bármelyik napon.

head(szeles[szeles$"102" > 0,c("store_nbr","date","102")])
##       store_nbr       date 102
## 36454         3 2013-03-08   1
## 36515         3 2013-03-29   1
## 36569         3 2013-05-06   1
## 36603         3 2013-05-17   1
## 36634         3 2013-06-22   1
## 36696         3 2013-07-30   1

Az érintett termékeket mindenhol kinullázom, ezzel ismét megoldottunk 140 sort. Egyszersmind törlöm a szóban forgó oszlopokat az adattáblából.

test.munka[test.munka$item_nbr %in% c("102","24","63","66"),"units"] <- 0
nrow(test.munka[is.na(test.munka$units),])
## [1] 26030
for (i in c("102","24","63","66")) {szeles[i] <- NULL}

3.3. Ünnepnapok kezelése

Az adatok szemrevételezése révén megállapítottam, hogy a Walmart-üzletek a jelek szerint az év minden napján nyitva tartanak (és van forgalmuk), kivéve december 25-én, tehát karácsonykor. Ezért a tesztadathalmazba 2013. december 25-ére mindenhol 0-t írok.

test.munka[test.munka$date == "2013-12-25", "units"] <- 0
nrow(test.munka[is.na(test.munka$units),])
## [1] 26024

Ennek sajnos nem sok haszna van, csak 6 sort sikerült kiküszöbölni, de érdemes volt megpróbálni.

3.4. Áruk bevezetésének, kivezetésének kezelése

Fent a 40., 93. és 99. termék kapcsán láttuk, hogy még ha egy termék kapható is egy üzletben, akkor is előfordulhat, hogy nem a teljes időszakban kapható, hanem annak csak egy részében, vagy még ha kapható is, egyáltalán nem vásárolnak belőle hosszú, több hónapig tartó időszakokban. Például a 40. terméknél pont 2013/2014 fordulóján van egy több hetes időszak, amikor 0 darab fogyott, és a 99. terméknél is látható egy hasonló időszak 2014 közepén (az adatok pontosabb vizsgálatából kiderül, hogy ez majdnem az egész június).

Ezek kezelésére a következő megoldást alkalmazom:

  1. Ha egy termékből egy adott boltban d dátumnál korábban nem adtak el egy darabot sem, akkor a d dátumnál korábbi napokra 0 darabot írok a tesztadatsorba. Ez csak akkor érdekes, ha d 2013. április után esik.
  2. Hasonlóan járok el a fordított esetben, tehát ha egy adott boltban d dátumnál később nem adtak el egy darabot sem egy árucikkből. Itt csak 2014. október 26. előtti dátumokat vizsgálom, mert odáig terjed a teszthalmaz.
  3. Az átmeneti készlethiányok (illetve csak szezonálisan forgalmazott árucikkek) kezelése gondot okoz. Az ugyanis, hogy két időpont között nem szerepel az eladások között egyetlen darab termék sem, egyrészt jelentheti azt, hogy ténylegesen nem történt eladás, másrészt azt is, hogy időjárási esemény(ek) miatt a kérdéses időszak nem szerepel a tanuló adathalmazban, és ugyan volt eladás, de nem szerepel a táblázatban. Ugyanakkor az is lényeges szempont, hogy valószínűleg sokkal több készlethiányos időszak van, mint ahány viharos időszak az egyes boltok területén. Ezért úgy járok el, hogy végigjárom boltonként a viharos időszakokat, és megnézem, hogy az előző és a következő napon vajon 0 szinten vannak-e az eladások. Ha igen, akkor kinullázom az adott viharos időszakot.

Természetesen mindhárom részmegoldás csupán heurisztikus jellegű. Az 1. és 2. eset viszonylag nem problematikus, ugyanis igen valószínűtlennek tűnik, hogy bármilyen termékből a látható bevezetése előtt, illetve kivezetése után kizárólag viharos napokon adtak volna el példányokat. Ennél nagyobb gondot okoz a 3. eset, ott bizonyára bizonyos mennyiségű hibával számolhatunk a vázolt egyszerű megoldás miatt (különösen olyan esetekben, ahol hosszú ideig folyamatosan tart egy viharos időszak, mint pl. a 9. állomásnál 2014. júniusban), azonban jobb megoldást nem látok.

for (b in c(1:34,36:45)) {
    s <- szeles[szeles$store_nbr==as.character(b),c(4:110)]
    s <- s[,colSums(s) != 0]
    s$date <- szeles[szeles$store_nbr==as.character(b),2]
    viharos <- test.munka[test.munka$store_nbr==as.character(b)
                          & test.munka$item_nbr==names(s[1]),]
    d <- viharos[1,1]
    elso <- d
    k <- 1
    while (k < nrow(viharos)) {
        while (k < nrow(viharos) & viharos[k+1,1] == viharos[k,1] + 1) {
            k <- k + 1
        }
        utolso <- viharos[k,1]
        for (cikk in names(s[1:(ncol(s)-1)])) {
            if (length(szeles[szeles$date == elso - 1 & szeles$store_nbr == b,cikk]) == 1 
                & length(szeles[szeles$date == utolso + 1 & szeles$store_nbr == b,cikk]) == 1)
                {
                if (szeles[szeles$date == elso - 1 
                         & szeles$store_nbr == b,cikk] == 0
                & szeles[szeles$date == utolso + 1 
                         & szeles$store_nbr == b,cikk] == 0) {
                test.munka[elso <= test.munka$date 
                           & test.munka$date <= utolso
                           & test.munka$store_nbr == b
                           & test.munka$item_nbr == cikk,"units"] <- 0
                }
            }
        }
        k <- k + 1
        elso <- viharos[k,1]
    }
    for (i in 1:(ncol(s)-1)) {
        sorok <- s[c(ncol(s),i)]
        cikk <- names(sorok[2])
        sorok <- sorok[sorok[2]>0,1]
        kezdet <- sorok[1]
        veg <- sorok[length(sorok)]
        if (kezdet > "2013-04-01") {
            test.munka[test.munka$date < kezdet
                       & test.munka$item_nbr == cikk,
                       "units"] <- 0            
        }
        if (veg < "2014-10-26") {
            test.munka[test.munka$date > veg
                       & test.munka$item_nbr == cikk,
                       "units"] <- 0            
        }
    }
    
}
nrow(test.munka[is.na(test.munka$units),])
## [1] 14316

Láthatóan ezekkel a lépésekkel csaknem megfeleztük a megbecsülendő sorok számát.

3.5. A hét napjainak és a hónapoknak a hatása

Számíthatunk arra, hogy a lakosság időbeosztása miatt a hétnek nem minden napján van az üzletekben azonos forgalom. Vizsgáljuk meg, hogy ez hogyan nyilvánul meg az összes vizsgált üzletben.

A tanuló adathalmazhoz hozzáadok egy oszlopot, amely a hét napját tartalmazza a dátum alapján. Az a változóba generálok egy táblázatot, amely azt tartalmazza, hogy az egyes árucikkekből mennyit adtak el a hét adott napján.

szeles$nap <- weekdays(szeles$date)
a <- aggregate(szeles[4:110],list(szeles$nap),sum)[c(2,3,5,1,4,6,7),]
names(a)[1] <- "nap"
a[1:12]
##         nap   1  10 100 101 103 104 105 106 107 108 109
## 2     hétfő 121  92  92  69  48  62  57  10  80 288 681
## 3      kedd 106 107  88  50  27  57  57  10  86 222 570
## 5    szerda 111 108  71  46  33  95  56  12  42 208 453
## 1 csütörtök  75  86  78  55  33  75  69  17  76 271 577
## 4    péntek  86 113  79  57  42 121 109  36  85 347 586
## 6   szombat 113 125  71  96  34 166 124  40  84 337 651
## 7  vasárnap 166 114  73  77  24  80  88  14  81 372 690

Ezután összeadom a sorokban szereplő számokat, és a következő eredményt kapom.

a$osszeg <- apply(a[2:108],1,sum)
a$szorzo <- a$osszeg / a$osszeg[4]
a[c("nap","osszeg","szorzo")]
##         nap osszeg   szorzo
## 2     hétfő 664290 1.202648
## 3      kedd 582677 1.054894
## 5    szerda 555945 1.006498
## 1 csütörtök 552356 1.000000
## 4    péntek 611608 1.107271
## 6   szombat 731490 1.324309
## 7  vasárnap 820325 1.485138

Láthatjuk, hogy hétfőn és hétvégén jelentősen több fogy az árukból összességében, mint csütörtökön. Ezért a modellezés előtt az összes eladási adatot beosztom a táblázatban szereplő szorzóval a hét napjának a függvényében, azaz az eladási adatokat normálom a csütörtöki szintre. A regresszió eredményeit pedig ugyanezzel a számmal fogom majd felszorozni.

Kérdés, hogy a hónapoktól is hasonlóan függ-e az eredmény, például a Walmartok környéke kiürül-e a nyári kánikula idején, vagy épp fordítva, mindenki a boltban hűsöl, és így fellendül a forgalom.

o <- szeles
o$honap <- months(o$date)
h <- aggregate(o[4:110],list(o$honap),sum)[c(5,4,9,1,8,7,6,2,12,11,10,3),]
h[1:12]
##       Group.1   1  10 100 101 103 104 105 106 107 108 109
## 5      január   0  85  80  31  45  14  13   7  76 117 457
## 4     február  31  74  42  37  17  11  45   1  52 160 331
## 9     március 203  99  71  64  24   6  14   1  70 183 388
## 1     április 130  64  43  49  19  20  14  13  54 225 428
## 8       május 109 103  30  24  10  85  53  13  38 224 268
## 7      június 128  58  50  13  23  77  56  31  18 207 254
## 6      július 145  77  22  24  16 129  88  23  17 105 194
## 2   augusztus  32  52  34  43  20 116  97  10  36 208 331
## 12 szeptember   0  26  25  30  12  60  64  26  23  60 254
## 11    október   0  41  55  39  21  57  44   7  50 315 459
## 10   november   0  29  48  35  25  47  36   5  46  93 390
## 3    december   0  37  52  61   9  34  36   2  54 148 454

Az már ebből a rövid részletből is kiderül, és nem is ér minket váratlanul, hogy az egyes árucikkek forgalmában jelentős szezonális ingadozások figyelhetőek meg, és ezek árucikkenként különbözőek, így kérdéses, hogy mennyire lehet célravezető egy olyanfajta normálás, amilyet a hét napjaival kapcsolatban hajtottam végre.

h$osszeg <- apply(h[2:108],1,sum)
h$szorzo <- h$osszeg / h$osszeg[11]
h[c("Group.1","osszeg","szorzo")]
##       Group.1 osszeg   szorzo
## 5      január 481121 1.702348
## 4     február 416366 1.473226
## 9     március 432796 1.531360
## 1     április 377516 1.335763
## 8       május 354982 1.256031
## 7      június 371423 1.314204
## 6      július 373337 1.320976
## 2   augusztus 405861 1.436056
## 12 szeptember 354522 1.254403
## 11    október 365986 1.294966
## 10   november 282622 1.000000
## 3    december 302159 1.069128

Érdekes módon az adatok szerint novemberben és decemberben kiugróan alacsony a forgalom, januárban kiugróan magas, míg az év maradék részében nagyjából egyenletes. Nyilvánvaló, hogy ugyanazt a hét napjainál alkalmazott normalizálás itt nem működnek, például a fenti táblázatban az 1. és a 104. termék teljesen szembe megy ezekkel a tendenciákkal. Meg lehetne próbálni termékenként normalizálni, azonban az azért nem látszik járhatónak, mert a hónapok (a hét napjaival ellentétben) nem függetlenek az időjárási adatoktól. Így jobb híján annyit tudok tenni, hogy kinullázom azokat a sorokat a teszthalmazban, amelyek olyan hónapban vannak, amelyekben a táblázat szerint egyetlen példány sem kelt el az adott termékből (tehát pl. az 1. terméknél kinullázok mindent januárban és szeptembertől decemberig).

test.munka$ho <- months(test.munka$date)
for (b in 2:108) {
    sorok <- h[b]
    cikk <- names(sorok)
    for (honap in 1:12) {
        if (sorok[honap,1] == 0) {
            test.munka[test.munka$ho == h[honap,1]
                 & test.munka$item_nbr == cikk,
                 "units"] <- 0            
        }
    }
}
nrow(test.munka[is.na(test.munka$units),])
## [1] 14310

Ezzel minimálisan csökkent a becsülendő adatok száma.

3.6. Regressziós modellek

A kinullázható adatoknak ezzel a végére értünk. A továbbiakban valamilyen statisztikai módszer alkalmazásával kell megbecsülnünk az eladási számokat. Mivel nem osztályozási feladattal van dolgunk, hanem számértékeket keresünk, különböző módszerek (klaszterezés, neurális hálók, SVM-ek) nem alkalmazhatóak ésszerűen. A legésszerűbb az adatokat lineáris regresszióval kiszámítani, mégpedig boltonként és termékenként.

A regressziós modellezés előtt végrehajtom a fent jelzett módon az eladások hét napja szerinti normálását a teljes tanító adathalmazban.

a2 <- a[c("nap","szorzo")]
szeles <- merge(szeles,a2,by="nap")
#### A merge megváltoztatja a sorok sorrendjét, ezért
#### rendezem és sorszámozom a sorokat
szeles <- szeles[order(szeles$store_nbr,szeles$date),]
row.names(szeles) <- 1:nrow(szeles)

for (i in 5:111) {szeles[i] <- szeles[i] / szeles$szorzo}

Ezt követően végigjárom a teszthalmazt, és ha találok benne olyan sort, amelyhez még nem határoztam meg értéket, akkor az adott üzlethez és árucikkhez létrehozok egy lineáris modellt, majd a kapott modellel kiszámítom az adott üzlethez és árucikkhez tartozó összes hiányzó adatot. Ha a modellben szereplő valamely változó értéke hiányzik az előrejelzendő sorhoz tartozó időjárási adatok közül, az összes előrejelzett érték átlagát adom meg eladási számként.

for (sor in 1:nrow(test.munka)) {
    if(is.na(test.munka[sor,"units"])) {
        b <- as.character(test.munka[sor,"store_nbr"])
        cikk <- as.character(test.munka[sor,"item_nbr"])
        tanit <- szeles[szeles$store_nbr == b,c(112:114,116:118)]
        tanit$y <- szeles[szeles$store_nbr == b,cikk]
        tanit$date <- szeles[szeles$store_nbr == b,"date"]
        
        ### A modellek illesztésekor figyelmen kívül hagyom azokat
        ### az időszakokat, amelyek a forgalmazás kezdete előtt vagy
        ### vége után vannak. Az átmeneti készlethiánnyal és a
        ### szezonalitással itt az egyszerűség kedvéért nem foglalkozom.
        kezdet <- tanit[tanit$y > 0,"date"][1]
        veg <- tanit[tanit$y > 0,"date"][length(tanit[tanit$y > 0,"y"])]
        tanit <- tanit[tanit$date >= kezdet,]
        tanit <- tanit[tanit$date <= veg,]
        
        ### Az illesztéskor figyelmen kívül hagyom azokat a sorokat,
        ### ahol az előrejelző változók bármelyike NA.
        tanit <- tanit[!is.na(rowSums(tanit[1:3,5:6])),]
        
        ### A hóesést csak akkor kezelem előrejelző változóként,
        ### ha az adott állomásnál esett egyáltalán hó, és nincsenek
        ### a hóesés oszlopban NA értékek.
        if (!is.na(sum(tanit$snowfall)) & sum(tanit$snowfall) > 0) {
            modell <- lm(y ~ tmax + tmin + tavg + snowfall 
                    + preciptotal + avgspeed, tanit)
        }
        else {
            modell <- lm(y ~ tmax + tmin + tavg + preciptotal + avgspeed, tanit)
        }
        alkalmaz <- test.munka[test.munka$store_nbr == b 
            & test.munka$item_nbr == cikk & is.na(test.munka$units),]
        alkalmaz <- merge(alkalmaz,key)
        alkalmaz <- merge(alkalmaz,weather,by=c("station_nbr","date"))
        alkalmaz$y <- predict(modell, alkalmaz)
        alkalmaz[is.na(alkalmaz$y),"y"] <- mean(alkalmaz$y,na.rm=TRUE)
        test.munka[test.munka$store_nbr == b 
            & test.munka$item_nbr == cikk
            & is.na(test.munka$units),"units"] <- alkalmaz$y
    }
}
nrow(test.munka[is.na(test.munka$units),])
## [1] 0

Láthatjuk, hogy sikeresen eltűnt az összes NA érték.

Végül a kapott értékeket felszorzom a normalizáló szorzó reciprokával, kerekítek, és a regresszió által előrejelzett negatív értékeket átírom 0-ra.

test.munka$nap <- weekdays(test.munka$date)
test.munka <- merge(test.munka,a2,by="nap")
test.munka$units <- test.munka$units * test.munka$szorzo
test.munka$units <- round(test.munka$units)
test.munka[test.munka$units < 0,"units"] <- 0

Törlöm a feldolgozás céljából felvett oszlopokat a test.munka táblázatból, majd egyesítem a test.munka és a test táblázatokat, és rendezem a sorokat az eredeti test.csv-nek megfelelő sorrendben.

for (i in c("nap","ho","szorzo")) {test.munka[i] <- NULL}
test <- rbind(test,test.munka)
test <- test[order(test$date,as.numeric(as.character(test$store_nbr)),as.numeric(as.character(test$item_nbr))),]

A versenyre az eredményeket egy olyan táblázatban kell leadni, amelynek két oszlopa van: egy store_nbr++item_nbr++date formátumú id mező és a units. Legenerálom az id mezőt, majd mentem az előrejelzéseket.

test$id <- paste(test$store_nbr,test$item_nbr,test$date, sep="_")
head(test)
##               date store_nbr item_nbr units             id
## 1       2013-04-01         2        1     0 2_1_2013-04-01
## 2       2013-04-01         2        2     0 2_2_2013-04-01
## 3       2013-04-01         2        3     0 2_3_2013-04-01
## 4       2013-04-01         2        4     0 2_4_2013-04-01
## 5045851 2013-04-01         2        5    62 2_5_2013-04-01
## 6       2013-04-01         2        6     0 2_6_2013-04-01
write.csv(test[c(5,4)],file="f:/downloads/submit.csv",row.names=FALSE,quote=FALSE)

4. Kiértékelés

Az eredményeimet a Kaggle-verseny automatikus értékelőjével értékeltettem ki (mivel a tényleges adatsort, amelyhez közelíteni kellett, nem tette közzé a Walmart). Az értékelés mutatószámaként itt az átlagos logaritmikus négyzetes hiba gyökét használták, lásd a Kaggle-verseny megfelelő oldalán.

Az elemzésem eredménye 128. lett, tehát a középmezőnyben végzett volna, mivel kb. 360 versenyző ért el érdemi, a teljesen kinullázott táblázat benyújtásához képest jobb eredményt. Sajnos a verseny határideje után három héttel nyújtottam be, ezért a helyezés nem érvényes. Ez ahhoz képest, hogy csak egy próbálkozásból jött ki, nem rossz eredmény, konkrétan a legjobb ilyen (a következő egyszer próbálkozó 179. lett, míg az első tíz helyezett átlagosan 68-szor próbálkozott). Ugyanakkor a Walmart ezzel az eredménnyel vélhetően nem kapkodott volna utánam az adatelemzői pozícióra.

pngfile <- "f:/downloads/helyezes.png"
pic <- readPNG(pngfile)
grid::grid.raster(pic)

Az eredmények javítása érdekében célszerű lehet megpróbálkozni nem lineáris regressziós modellekkel, bár a tanító adathalmaz szemrevételezése alapján úgy tűnik, hogy ez sem sokat segített volna, mivel nagyon sok a zaj az adatokban. Ennél ígéretesebbnek tűnik annak vizsgálata, hogy az azonos kódszámú árucikkek eladásai a különböző üzletekben hasonlóan alakultak-e. Ha kiderülne, hogy igen, akkor azon üzletek eladási adatai alapján, ahol éppen nincs vihar, valószínűleg pontosabban közelíteni lehetne az eladások alakulásához a vihar által érintett üzletekben.