Karetní transakce na čerpacích stanicích

Úvod

Cílem této práce je zpracovat data, kde banka sbírala údaje o tom, jak lidé platili kartou na pumpách především v České Republice, ale i v zahraničí a to v letech 2017-2018.

Popis datasetu

Jde o csv soubor, kde každy řádek zaznamenává jednu transakci. Dataset obsahuje následující sloupce.

clid: identifikace klienta
posid: identifikace čerpací stanice
day: datum platby ve formátu Y-m-d
hour: hodina platby (0-23)
amt: částka platby v CZK
fav: zda je tato stanice označena jako oblíbená klienta
fav_main: zda je tato stanice označena jako nejoblíbenější klienta
cl_gps_lat, cl_gps_lon: přibližné GPS souřadnice bydliště klienta
country: země platby
city: město platby
name: název čerpací stanice
clst: číslo clusteru, do kterého stanice spadá na základě chování klientů (0 = nezařazeno do žádného clusteru)
pos_gps_lat, pos_gps_lon: GPS souřadnice místa platby (platí pouze pro stanice v České republice)
dist: vzdálenost v km mezi bydlištěm klienta a místem platby

Csv soubor má 1853551 zázamů. Jeho velikost je 226 MB.

Cíle

Důvodem této analýzy je především získání zápočtu a dále snaha co nejlépe porozumět lidskému chování při placení na čerpacích stanicích na základě těchto dat. Tato analýza může mít užitek pro čerpací stanice, které mohou na základě dat optimalizovat své chování a zjistit, jak generovat větší zisky nebo jaká očekávání mohou mít pro dané období. Navíc i pro banky mohou být tato data cenná, neboť mohou zjistit, jakým zákazníkům je rozumné nabízet služby spojené například s cestováním, pojistkou na auto a jiné s tímto spojené.

Analýza vlastností datasetu

V této části více prozkoumám dataset s cílem získat co nejvíce vhledu do toho jak v něm data vypadají a jaké jsou jeho nedokonalosti (chybející / špatné hodnoty).

Struktura datasetu

Takto vypadají řádky - zde je náhodný sample 5ti z nich:

Code

import pandas as pd

file_path = 'benz-tr.csv'
df = pd.read_csv(file_path)

df.sample(5)

	clid	posid	day	hour	amt	fav	fav_main	cl_gps_lat	cl_gps_lon	country	city	name	clst	pos_gps_lat	pos_gps_lon	dist
1490060	19892481	8059	2017-12-15	16.0	1128.98	1	1	49.677821	13.027792	cze	stribro	mol stribro 709	1.0	49.751507	12.998950	8.451807
544606	10330931	5693	2017-08-07	20.0	37.00	1	1	50.634210	15.489332	cze	jilemnice	gbo hrabacov	8.0	50.643860	15.197062	20.638822
1561060	18214895	8277	2017-07-11	10.0	97.00	1	1	49.978352	16.569073	cze	trutnov	benzina cs 0414	1.0	50.572090	15.893961	81.611036
250934	9864640	4908	2017-11-20	14.0	310.00	0	0	50.303806	12.561612	cze	chodov	benzina cs 0108	2.0	50.239868	12.754909	15.468209
957318	18569452	6904	2017-12-01	14.0	159.80	0	0	49.626580	18.073860	cze	opava jaktar	mol opava 008	7.0	49.952827	17.869177	39.139811

Podívejme se kolik různých klientů bylo sledováno a na kolika různách pumpách platili.

Code

unique_clients = df['clid'].nunique()
unique_posids = df['posid'].nunique()

unique_df = pd.DataFrame({
    'Počet klientů': [unique_clients],
    'Počet pump': [unique_posids]
})

unique_df

	Počet klientů	Počet pump
0	121834	23856

Dataset obsahuje duplicity

Code

duplicate_rows = df[df.duplicated(keep=False)]
print("Duplicitních řádků je: ", duplicate_rows.shape[0])
# non_czech_dup = duplicate_rows[duplicate_rows["country"] != "cze"]
# non_czech_dup.shape[0]

Duplicitních řádků je:  1723

Bohužel nevím jak dataset vznikal a jak lépe interpretovat existenci těchto duplikátů. Je možné, že i některé nebudou vyloženě chybné duplikáty. Některé mohli vzniknout tak, že klient chvíli po první transakci na dané benzínce provedl další za stejnou sumu (například si klient vzpomněl, že i jeho žena chtěla koupit bagetu). Ovšem většina duplicitních transakcí obsahuje vyšší částky a tak určitě duplikáty z důvodu opakovaného nákupu nebudou tak časté a půjde spíše opravdu o dvojitý záznam zvniknuvší nějakou chybou.

Code

duplicate_rows.sample(5)

	clid	posid	day	hour	amt	fav	fav_main	cl_gps_lat	cl_gps_lon	country	city	name	clst	pos_gps_lat	pos_gps_lon	dist
1735706	16370954	8786	2017-06-24	8.0	19.00	0	0	49.650314	18.022462	cze	zebrak	cerpaci stanice kapr	7.0	49.869896	13.894603	297.472603
1044559	15214318	7060	2018-01-12	17.0	86.00	1	1	49.829670	18.253827	cze	ostrava	omv 2119	11.0	49.804031	18.287262	3.725882
1801943	13379	28204	2017-09-09	12.0	44.01	0	0	49.837316	18.279926	pol	plock	orlen	NaN	NaN	NaN	NaN
1853114	16455167	42014	2018-01-09	15.0	1053.38	0	0	49.085408	14.719735	swe	vesteras	ingo vesteras stromle	NaN	NaN	NaN	NaN
1317998	18022577	7608	2017-11-16	6.0	84.00	0	0	50.050808	14.489218	cze	praha	omv 2129	13.0	50.014519	14.425369	6.089410

Nevím jestli do tohoto datasetu byly evidovány i transakce jež skončily chybou. Tato možnost mi dává smysl, protože značná část duplicit je ze zahraničí, kde spíše proběhne nějaký problém při platbě a je třeba ji opakovat. Lze si také povšimnout, že hodně duplicit má častku dělitelnou číslem 100. K tomuto se dostanu později.

Code

duplicate_rows_outside_cze = duplicate_rows[duplicate_rows["country"] != "cze"]
total_amount_of_rows_outside_cze = df[df["country"] != "cze"].shape[0]

duplicate_rows_divisible_by_100 = duplicate_rows[duplicate_rows["amt"] % 100 == 0]
total_amount_of_rows_divisible_by_100 = df[df["amt"] % 100 == 0].shape[0]

print(f"Duplicitních plateb, které proběhly mimo Českou Republiku je {duplicate_rows_outside_cze.shape[0] / duplicate_rows.shape[0] * 100:.2f}%. Ovšem mimo ČR celkově v datasetu proběhlo jen {total_amount_of_rows_outside_cze / df.shape[0] * 100:.2f}% plateb.")
print(f"Duplicitních plateb, které jsou dělitelné číslem 100 je {duplicate_rows_divisible_by_100.shape[0] / duplicate_rows.shape[0] * 100:.2f}%.")

Duplicitních plateb, které proběhly mimo Českou Republiku je 18.40%. Ovšem mimo ČR celkově v datasetu proběhlo jen 4.41% plateb.
Duplicitních plateb, které jsou dělitelné číslem 100 je 22.17%.

Chyby v datasetu

Bohužel jak tomu často bývá, i v tomto datasetu jsou chybějící údaje. Zde je tabulka v jakých sloupcích data chybí a kolik jich chybí:

Code

missing_values = df.isnull().sum()
total_rows = df.shape[0]
missing_percentage = (missing_values / total_rows) * 100

missing_data_analysis = pd.DataFrame({
    'Chybějících hodnot': missing_values,
    'Percento z celkového počtu': missing_percentage
})

missing_data_analysis = missing_data_analysis[missing_data_analysis['Chybějících hodnot'] > 0]
missing_data_analysis

	Chybějících hodnot	Percento z celkového počtu
hour	2631	0.141944
cl_gps_lat	14413	0.777589
cl_gps_lon	14413	0.777589
country	10	0.000540
city	135	0.007283
clst	254440	13.727165
pos_gps_lat	117187	6.322297
pos_gps_lon	117187	6.322297
dist	130214	7.025110

Vidíme, že cl_gps_lat a lon chybí, tak jak by jeden očekával, vždy společně a to samé platí pro pos_gps_lat a lon. (Byť jich chybí stejně, tak toto není nutně důkaz, ale lze to ověřit i přes pandas.)

Code

missing_cl_gps = df[df['cl_gps_lat'].isnull() & df['cl_gps_lon'].isnull() & df['cl_gps_lon'].isnull()]
all_missing_cl_gps = missing_cl_gps.shape[0] == df['cl_gps_lat'].isnull().sum()

missing_pos_gps = df[df['pos_gps_lat'].isnull() & df['pos_gps_lon'].isnull() & df['pos_gps_lon'].isnull()]
all_missing_pos_gps = missing_pos_gps.shape[0] == df['pos_gps_lat'].isnull().sum()

# all_missing_cl_gps, all_missing_pos_gps

Podívejme se na korelace, kdy společne chybí hodnoty. Jelikož mi přišlo i zajimavé sledovat, zdalipak země transakce byla mimo ČR, tak jsem přidělal dummy sloupec, aby se to zahrnulo také do korelace.

Code

df2 = df.copy()
df2['Is_not_cze'] = 0
df2.loc[df['country'] != 'cze', 'Is_not_cze'] = None
missing_value_correlation = df2.isnull().corr()

high_corr_pairs = missing_value_correlation.unstack().sort_values(kind="quicksort", ascending=False)
high_corr_pairs = high_corr_pairs[high_corr_pairs < 1]  # Remove self-correlation (which is always 1)

high_corr_pairs_df = high_corr_pairs.reset_index()
high_corr_pairs_df.columns = ['Sloupec 1', 'Sloupec 2', 'Míra korelace']
high_corr_pairs_df = high_corr_pairs_df.drop_duplicates(subset=['Míra korelace'])

high_corr_pairs_df[high_corr_pairs_df['Míra korelace'] > 0.1]

	Sloupec 1	Sloupec 2	Míra korelace
0	pos_gps_lon	dist	0.945095
4	pos_gps_lon	Is_not_cze	0.826345
8	Is_not_cze	dist	0.780975
10	Is_not_cze	clst	0.538180
12	clst	pos_gps_lon	0.461532
16	clst	dist	0.437353
18	cl_gps_lon	dist	0.322052

Je vidět, že ujetá vzádelost chybí, když nevíme lokalitu bydliště člověka nebo benzínky, což je logické, protože pak není jak jí spočítat. Zařazení do clusteru není provedeno opět když nevíme bydliště či pozici benzínky. Tyto data chybí převážně, když země transakce není Česko. Ovšem další chybějící data vyžadují podrobnější analýzu kdy chybí a proč.

Podívejme se na to kdy chybí hodiny.

Code

missing_hours_cze = df[
    (df['country'] == 'cze') & df['hour'].isnull()
].shape[0]

missing_hours_not_cze = df[
    (df['country'] != 'cze') & df['hour'].isnull()
].shape[0]

total_missing_hours = df['hour'].isnull().sum()

print(f"V ČR chybí {missing_hours_cze} hodinových údajů, mimo ČR chybí {missing_hours_not_cze} hodinových údajů. Celkem chybí {total_missing_hours} hodinových údajů.")
print("Zde je tabulka s deseti nahodnými samply, kdy chybí hodinový údaj:")
df[df['hour'].isnull()].sample(7)

V ČR chybí 2464 hodinových údajů, mimo ČR chybí 167 hodinových údajů. Celkem chybí 2631 hodinových údajů.
Zde je tabulka s deseti nahodnými samply, kdy chybí hodinový údaj:

	clid	posid	day	hour	amt	fav	fav_main	cl_gps_lat	cl_gps_lon	country	city	name	clst	pos_gps_lat	pos_gps_lon	dist
1298881	14395982	7558	2017-02-08	NaN	1655.80	0	0	50.000017	14.560070	cze	pruhonice	union road cerpaci st	16.0	50.000007	14.559255	0.058276
1833516	18014522	37680	2017-02-21	NaN	399.97	0	0	50.030059	13.929977	cze	krivoklat	cerpaci stanice kri	NaN	50.037201	13.877518	3.830136
453743	15844343	5462	2017-02-28	NaN	300.00	1	0	50.184745	15.835921	cze	hradec kralove	benzina cs 0446	16.0	50.202755	15.840381	2.027638
1837499	15231755	37756	2017-03-07	NaN	500.00	1	1	49.849272	18.001223	cze	ostrava	globus tankautomat	16.0	49.846191	18.148991	10.600599
1833967	14329874	37694	2017-02-22	NaN	500.00	1	1	50.064019	16.551233	cze	letohrad	samoobsluzna cerpaci	16.0	50.031953	16.507928	4.719604
1834690	18131802	37703	2017-02-16	NaN	500.00	1	1	50.557408	14.032168	cze	litomerice	cargonet octomat	16.0	50.539491	14.129419	7.154374
1845219	15166461	37827	2017-02-07	NaN	282.60	0	0	49.333379	16.478677	cze	pribram	ps cng pribram	NaN	49.679956	13.999663	183.090660

Nejedná se tedy o chyby, jak jsem předpokládal, že například platba bude v jiné zemi a bude to kvůli rozdílnému časovému pásmu, či že zkrátka jde o obyčejné nahodné chyby, ale souvisí to převážně s divnými transakcemi, jež jsou hezká čísla (dělitelná stem). K těmto transakcím se dostanu později.

Code

missing_hour_df = df[df['hour'].isnull()]

divisible_by_100 = missing_hour_df['amt'] % 100 == 0
divisible_by_100_count = divisible_by_100.sum()

divisible_by_100_df = pd.DataFrame({
    'Celkem chybějících záznamů hodin': [missing_hour_df.shape[0]],
    'Počet řádků, kde chybí hodina a platba je dělitelná stem': [divisible_by_100_count]
})

divisible_by_100_df

	Celkem chybějících záznamů hodin	Počet řádků, kde chybí hodina a platba je dělitelná stem
0	2631	1504

Code

print(f"Celkem tedy {divisible_by_100_count / missing_hour_df.shape[0] * 100:.2f}% záznamů má částku dělitelnou stem, když chybí hodina.")

Celkem tedy 57.16% záznamů má částku dělitelnou stem, když chybí hodina.

Údaj o městu chybí na jedné beznínce na Slovensku a v Řecku. V obou případech jde o Shellku.

Code

missing_cities_df = df[df['city'].isnull()]

if not missing_cities_df.empty:
    missing_cities_by_country = missing_cities_df['country'].value_counts()
    top_countries_missing_cities_df = pd.DataFrame({
        'Země': missing_cities_by_country.index,
        'Chybějících údajů s městem': missing_cities_by_country.values,
        "Pos_idčka pump": missing_cities_df["posid"].unique()
    })


# print(missing_cities_df[missing_cities_df["posid"] == 31271].shape[0])
top_countries_missing_cities_df

	Země	Chybějících údajů s městem	Pos_idčka pump
0	svk	132	19507
1	grc	3	31271

V datasetu chybí 10 záznamů o zemi. Zemi, ale lze určit z názvu města, kde se pumpa nachází - jedná se vždy o Kosovo (pouze částečně uznaný stát, možná proto tam není nic). Hodnoty tak doplním.

Code

missing_coutry = df[df["country"].isna()].groupby("city").count()["posid"].sort_values(ascending=False)
df.loc[df["city"].isna(), "city"] = "xk"
missing_coutry

city
prishtine    5
gjakove      2
besiana      1
peje         1
prizren      1
Name: posid, dtype: int64

Exploratorní analýza dat

V této části se ponořím hlouběji do analýzy dat, s tím že nyní budou i koukat na vzorce které vykazují a.j.

Zde je basic overview statistik sloupců:

Code

pd.set_option('display.float_format', lambda x: '%d' % x if x.is_integer() else '%.2f' % x)
distribution_stats = df.describe()
unique_entries = df.nunique()
distribution_stats.loc['unique'] = unique_entries
distribution_stats

	clid	posid	hour	amt	fav	fav_main	cl_gps_lat	cl_gps_lon	clst	pos_gps_lat	pos_gps_lon	dist
count	1853551	1853551	1850920	1853551	1853551	1853551	1839138	1839138	1599111	1736364	1736364	1723337
mean	15470252.61	7608.88	13.39	551.17	0.55	0.39	49.87	15.60	5.46	49.85	15.59	26.37
std	3044302.51	5242.39	4.71	517.25	0.50	0.49	0.50	1.59	4.81	0.50	1.58	49.37
min	8081	1	0	10	0	0	48.60	12.16	0	48.59	12.17	0
25%	14243273	5478	10	160	0	0	49.50	14.42	1	49.47	14.44	2.50
50%	15735517	6857	14	400.10	1	0	49.92	15.13	5	49.88	15.16	7.61
75%	18043012	7786	17	865.70	1	1	50.17	16.76	9	50.16	16.71	22.05
max	22410503	42750	23	38990	1	1	51.02	18.83	17	51.03	18.76	455.75
unique	121834	23856	24	118776	2	2	8545	8545	18	2140	2147	247147

Analýza plateb

Nyní se zaměřme na platby. Zde je graf, jenž ukazuje distribuci částek.

Code

import seaborn as sns
import matplotlib.pyplot as plt

sns.set_style("whitegrid")

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(18, 9))

sns.histplot(df[df["amt"] < 2500]["amt"], kde=True, bins=100, ax=axs[0])
axs[0].set_title('Distribuce plateb')
axs[0].set_xlabel('Výška platby')
axs[0].set_ylabel('Frekvence')

sns.ecdfplot(data=df[df["amt"] < 2500], x="amt", ax=axs[1])
axs[1].set_title('Kumulativní distribuce plateb')
axs[1].set_xlabel('Výška platby')
axs[1].set_ylabel('Kumulativní proporce')

plt.tight_layout()
plt.show()

Je vidět, že většina lidí na beznínce nechá okolo 500ti korun a naprostá většina se vejde do dvou tisíc. Opět lze tady pozorovat jak moc jsou časté platby jako je 300, 500, 1000 kč.

Code

top5_common_amt_percentage = df['amt'].value_counts(normalize=True).head(5)*100
top5_common_amt_df = top5_common_amt_percentage.reset_index()
top5_common_amt_df.columns = ['Platba (CZK)', 'Procento z celkových plateb']
top5_common_amt_df

	Platba (CZK)	Procento z celkových plateb
0	500	2.12
1	300	1.32
2	86	1.27
3	1000	1.10
4	200	1.07

Code

payments_divisible_by_100 = df[df['amt'] % 100 == 0]


payments_by_brand = payments_divisible_by_100.groupby('name')['amt'].count().reset_index().sort_values(by='amt', ascending=False).head(10)
payments_by_brand.columns = ['Jméno benzínky', 'Počet plateb dělitelných 100']

payments_by_brand

	Jméno benzínky	Počet plateb dělitelných 100
790	cs eurooil	17489
427	benzina cs 0616	2413
1235	globus tankautomat	1949
1134	csphm tank ono	1878
228	benzina cs 0279	1759
1923	robin oil	1664
326	benzina cs 0446	1570
2199	unixan cs unicorn	1541
2198	union road cerpaci st	1450
450	benzina cs 0644	1437

Bohužel přimo z dat nejde vyčíst důvod k tomu proč je tolik plateb takto “hezkých”. Nemyslím si, že by se jednalo o dýška, jelikož při platbě kartou na beznínce jej opravdu minimum lidí bude dávat. Je možné (tohle považuji jako za dost pravděpodobné), že jde o platby u samoobslužných pump, kde člověk musí prvně dát kartu a limit do kolika chce tankovat. Tomu i nasvědčuje to, že eurooil poskytuje tyto automaty, globus také atd. Bohužel toto ovšem přímo nevysvětluje proč čas chybí často u těchto hodnot - možná je to tím, že platba je zadaná předem, vytvoří se blokace na kartě a strhne se platba až později po natankování - tam můžou být prodlevy a tak se neviduje správně hodina. To by i částečně vysvětlovalo, proč duplicity mají také často tyto částky - často se to pokazí a také se to může zaevidovat 2x (zaeviduje se blokace a i ztrhnutí).

Podívejme se nyní na to jak se platí v ruzných ročních obdobích.

Code

df['day'] = pd.to_datetime(df['day'])

def get_season(month):
    if month in [12, 1, 2]:
        return 'Zima'
    elif month in [3, 4, 5]:
        return 'Jaro'
    elif month in [6, 7, 8]:
        return 'Podzim'
    else:
        return 'Léto'

df['season'] = df['day'].dt.month.apply(get_season)

season_avg_amount = df.groupby('season')['amt'].mean().sort_values(ascending=False)

season_total_amount = df.groupby('season')['amt'].count().sort_values(ascending=False)

colors1 = ['cyan', 'limegreen', 'red', 'orange']
colors2 = ["orange", "red", "limegreen", "cyan"]

plt.subplot(1, 2, 1)
season_avg_amount.plot(kind='bar', figsize=(18, 9), color=colors1, edgecolor='black', alpha=0.7)
plt.title('Průmerná výše platby podle ročního období', fontsize=16)
plt.xlabel('Roční období', fontsize=14)
plt.ylabel('Průměrná výše platby (CZK)', fontsize=14)
plt.grid(axis='y', linestyle='--', linewidth=0.5)
plt.xticks(rotation=45)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  

plt.subplot(1, 2, 2)
season_total_amount.plot(kind='bar', figsize=(18, 9), color=colors2, edgecolor='black', alpha=0.7)
plt.title('Počet plateb za období', fontsize=16)
plt.xlabel('Roční období', fontsize = 14)
plt.ylabel('Počet plateb', fontsize=14)
plt.grid(axis='y', linestyle='--', linewidth=0.5)
plt.xticks(rotation=45)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  

plt.show()

V prvním grafu vidíme, že nejvyšší průměrné platby jsou dělány v Zimě. Důvodem může být to, že auta spotřebovávají více benzínu, kvůli nižší efektivitě motoru za nízkých teplot nebo to, že kvůli zimě tankují klienti raději méně často (také kvůli Vánocům a silvestru) a pak tedy i za více.

Druhý graf není tolik vypovídající, protože data jsou sbírána od února do konce února dalšího roku a spíše odráží to, že banka nabírala nové klienty a tak roste i celkový počet evidonavaných plateb.

Code

unique_clients_upto_date = df.groupby('day')['clid'].apply(lambda x: x.unique())

unique_clients_set = set()
cumulative_counts = []

for date, clients in unique_clients_upto_date.items():
    unique_clients_set.update(clients)
    cumulative_counts.append(len(unique_clients_set))

cumulative_unique_clients_series = pd.Series(cumulative_counts, index=unique_clients_upto_date.index)

plt.figure(figsize=(18, 9))
cumulative_unique_clients_series.plot()
plt.title('Kumulativní počet unikátních klientů co provedli tržbu do daného data', fontsize=16)
plt.xlabel('Datum', fontsize=14)
plt.ylabel('Kumulativní počet unikátních klientů', fontsize=14)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=12)  
plt.grid(True)
plt.show()

Code

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(20, 8))

avg_payment_by_hour = df.groupby('hour')['amt'].mean().reset_index()
sns.barplot(data=avg_payment_by_hour, x='hour', y='amt', ax=axs[0], palette="viridis")
axs[0].set_title('Průměrná výše platby podle hodiny', fontsize=15)
axs[0].set_xlabel('Hodina', fontsize=14)
axs[0].set_ylabel('Průměrná výše platby', fontsize=14)
axs[0].set_xticklabels(axs[0].get_xticklabels(), rotation=45, fontsize=14)



count_payments_by_hour = df['hour'].value_counts().sort_index().reset_index()
count_payments_by_hour.columns = ['hour', 'count']
sns.barplot(data=count_payments_by_hour, x='hour', y='count', ax=axs[1], palette="magma")
axs[1].set_title('Počet plateb podle hodiny', fontsize=15)
axs[1].set_xlabel('Hodina', fontsize=14)
axs[1].set_ylabel('Počet plateb', fontsize=14)
axs[1].set_xticklabels(axs[1].get_xticklabels(), rotation=45, fontsize=14)

plt.tight_layout()
plt.show()

Výše plateb v pozdních večerních hodinách je na beznínkách nižší - nejspíše je lidé užívají hodně jako “večerky”, ale jak ukazuje druhý graf, tak se jedná pouze o malou část lidí co pozdě večer takto nakupují. Peak lidí můžou beznínky očekávat okolo 14-18 hodiny, kde i výška plateb je nejvyšší.

Code

total_amt_by_hour = df.groupby('hour')['amt'].sum().reset_index()

plt.figure(figsize=(18, 8))
sns.barplot(data=total_amt_by_hour, x='hour', y='amt', palette="coolwarm")
plt.title('Celkové držby benzínek dle hodiny', fontsize=16)
plt.xlabel('Hodina', fontsize=14)
plt.ylabel('Celková suma', fontsize=14)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  
plt.show()

Je důležité pro benzínky si spočítat, zdalipak se jim vyplatí operovat mezi 23-05 hodinou, kdy provoz potenicálně může být drazší než kolik vydělají.

Na nejvíce zákazníků se pumpy musí připravit v pátek, kdy lidé asi tankují na cesty např. na chatu kam jedou na víkend atd.

Code

df['day_of_week'] = df['day'].dt.day_name()

avg_payment_by_weekday = df.groupby('day_of_week')['amt'].mean()
total_payment_by_weekday = df.groupby('day_of_week')['amt'].sum()

order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
avg_payment_by_weekday = avg_payment_by_weekday.reindex(order)
total_payment_by_weekday = total_payment_by_weekday.reindex(order)

day_name_mapping = {
    'Monday': 'Pondělí',
    'Tuesday': 'Úterý',
    'Wednesday': 'Středa',
    'Thursday': 'Čtvrtek',
    'Friday': 'Pátek',
    'Saturday': 'Sobota',
    'Sunday': 'Neděle'
}

total_payment_by_weekday.index = total_payment_by_weekday.index.map(day_name_mapping)

plt.figure(figsize=(16, 8))

sns.barplot(x=total_payment_by_weekday.index, y=total_payment_by_weekday.values)
plt.title('Celková částka plateb podle dne v týdnu', fontsize=15)
plt.xlabel('Den v týdnu', fontsize=13)
plt.ylabel('Celková částka plateb', fontsize=13)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

plt.show()

Code

df['day_of_month'] = df['day'].dt.day

transactions_by_day = df.groupby('day_of_month').size()
avg_payment_by_day = df.groupby('day_of_month')['amt'].mean()

plt.figure(figsize=(18, 7))

plt.subplot(1, 2, 1)
transactions_by_day.plot(kind='bar', color='skyblue')
plt.title('Počet transakcí podle dne v měsíci', fontsize=15)
plt.xlabel('Den v měsíci', fontsize=13)
plt.ylabel('Počet transakcí', fontsize=13)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

plt.subplot(1, 2, 2)
avg_payment_by_day.plot(kind='bar', color='lightgreen')
plt.title('Průměrná výše platby podle dne v měsíci', fontsize=15)
plt.xlabel('Den v měsíci', fontsize=13)
plt.ylabel('Průměrná výše platby', fontsize=13)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
plt.tight_layout()
plt.show()

Analýza klientů

Zde je korelační matice klientů - z ní se dá zjistit jak jedno chovnání může ovlivnit jiné.

Code

df_dummies = df.copy()
df_dummies['payment_outside_cze'] = df['country'] != 'cze'
df_dummies['payment_inside_cze'] = df['country'] == 'cze'
df_dummies['payment_divisible_by_100'] = df['amt'] % 100 == 0

grouped_df_by_clients = df_dummies.groupby('clid').agg(
    total_spent=pd.NamedAgg(column='amt', aggfunc='sum'),
    avg_spent=pd.NamedAgg(column='amt', aggfunc='mean'),
    unique_stations=pd.NamedAgg(column='posid', aggfunc='nunique'),
    fav_stations=pd.NamedAgg(column='fav', aggfunc='mean'),
    fav_main_stations=pd.NamedAgg(column='fav_main', aggfunc='mean'),
    avg_distance=pd.NamedAgg(column='dist', aggfunc='mean'),
    count_payment_outside_cze=pd.NamedAgg(column='payment_outside_cze', aggfunc='sum'),
    count_payment_inside_cze=pd.NamedAgg(column='payment_inside_cze', aggfunc='sum'),
    count_payment_divisible_by_100=pd.NamedAgg(column='payment_divisible_by_100', aggfunc='sum'),
    amount_of_records=pd.NamedAgg(column='clid', aggfunc='count'),
    avg_payment_hour=pd.NamedAgg(column='hour', aggfunc='mean')
)

grouped_df_by_clients = grouped_df_by_clients.rename(columns={
    'total_spent': 'celkově utraceno',
    'avg_spent': 'průmerná útrata',
    'unique_stations': 'navštíveno unikátních stanic',
    'fav_stations': 'frekvence navštěvovaní oblíbených stanic',
    'fav_main_stations': 'frekvence navštěvovaní nejoblíbenější stanice',
    'avg_distance': 'průměrná vzdálenost',
    'count_payment_outside_cze': 'počet plateb mimo ČR',
    'count_payment_inside_cze': 'počet plateb v ČR',
    'count_payment_divisible_by_100': 'počet plateb delitelných 100',
    'amount_of_records': 'počet plateb udělaných klientem',
    'avg_payment_hour': 'průměrná hodina platby'
})

grouped_corr_matrix = grouped_df_by_clients.corr()

plt.figure(figsize=(16, 9))
sns.heatmap(grouped_corr_matrix, annot=True, cmap='coolwarm', fmt=".2f", annot_kws={"size": 10}, xticklabels=True, yticklabels=True)
plt.title('Korelační matice klientských dat', fontsize=20)
plt.xticks(fontsize=14, rotation=45, ha='right')
plt.yticks(fontsize=14)
plt.show()

Podívejme se nyní na to jak klienti různě platí a kolik navštívili různých stanic:

Code

grouped_df = df.groupby('clid').agg(
    number_of_payments=pd.NamedAgg(column='amt', aggfunc='size'),
    number_of_distinct_stations=pd.NamedAgg(column='posid', aggfunc='nunique'),
    average_amount=pd.NamedAgg(column='amt', aggfunc='mean'),
    average_hour=pd.NamedAgg(column='hour', aggfunc='mean')
)

plt.figure(figsize=(18, 30))

plt.subplot(4, 1, 1)
sns.histplot(grouped_df[grouped_df['number_of_payments'] < 50]['number_of_payments'], kde=True, color='blue', bins=30)  # Adjust the number of bins
plt.title('Distribuce počtu plateb na klienta', fontsize=20)
plt.xlabel('Počet plateb', fontsize=20)
plt.ylabel('Počet klientů', fontsize=20)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  

plt.subplot(4, 1, 2)
sns.histplot(grouped_df[grouped_df['number_of_distinct_stations'] < 30]['number_of_distinct_stations'], kde=True, color='green', bins=30)  # Adjust the number of bins
plt.title('Distribuce počtu navšívených různých stanic na klienta', fontsize=20)
plt.xlabel('Počet stanic', fontsize=20)
plt.ylabel('Počet klientů', fontsize=20)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  

plt.subplot(4, 1, 3)
sns.histplot(grouped_df[grouped_df['average_amount'] < 2000]['average_amount'], kde=True, color='red', bins=100)  # Adjust the number of bins
plt.title('Distribuce průměrné výše platby na klienta', fontsize=20)
plt.xlabel('Průměrná platba', fontsize=20)
plt.ylabel('Počet klientů', fontsize=20)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)

plt.subplot(4, 1, 4)
sns.histplot(grouped_df['average_hour'], kde=True, color='brown', bins=24)  # Adjust the number of bins
plt.title('Distribuce průměrné hodiny platby', fontsize=20)
plt.xlabel('Průměrný čas', fontsize=20)
plt.ylabel('Počet klientů', fontsize=20)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14) 

plt.tight_layout()
plt.show()

K mému překvapení se opět ukazuje, že nejoblíbenější čas na tankování je ve 13-14 hodin, v čase kdy bych očekával, že většina lidí pracuje, je možné, že se jedná o řidiče kamionů, či jiné pracující lidi, jejichž práce zahrnuje ježdení autem a kteří nyní mají touto dobou pauzu na oběd, a tak i k tomu tankují.

Podívejme s v jakých městech klienti nejčastěji tankují a tedy asi také i bydlí

Code

from wordcloud import WordCloud
df['city'] = df['city'].str.lower().str.strip()
text = " ".join(city for city in df.city.dropna() if (" nad " not in city and " u " not in city and " pod " not in city))

wordcloud = WordCloud(collocations=False, width=800, height=400, background_color='white').generate(text)
plt.figure(figsize=(14, 7))
plt.imshow(wordcloud, interpolation='bilinear')
plt.axis("off")
plt.show()

Cestují zakazníci v ruzných měsících jinak? To lze usoudit ze vzdálenosti mezi bydlištěm a místem platby. Nepřekvapivě peaky jsou v červenci a v srpnu, kdy rodiny často autem jedou do zahraničí na dovolenou. V lednu cestují lidé nejméně.

Code

df['month'] = pd.to_datetime(df['day']).dt.month

avg_distance_per_month_short = df.groupby('month')['dist'].mean()

plt.figure(figsize=(16, 8))
avg_distance_per_month_short.plot(kind='bar', color='orange', edgecolor='black', alpha=0.7)
plt.title('Průměrná vzdálenost dle měsíce', fontsize=20)
plt.xlabel('Měsíc', fontsize=16)
plt.ylabel('Průměrná vzdálenost (km)',fontsize=16)
plt.xticks(rotation=0)
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.show()

Analýza benzínek

V této části se obdovně jako v předchozí zaměřím na vzorce, které pro tentokrát vykazují benzínové pumpy.

Code

df_dummies = df.copy()
df_dummies['payment_outside_cze'] = df['country'] != 'cze'
df_dummies['payment_inside_cze'] = df['country'] == 'cze'
df_dummies['payment_divisible_by_100'] = df['amt'] % 100 == 0

grouped_df_by_stations = df_dummies.groupby('posid').agg(
    total_amount=pd.NamedAgg(column='amt', aggfunc='sum'),
    average_transaction=pd.NamedAgg(column='amt', aggfunc='mean'),
    number_of_clients=pd.NamedAgg(column='clid', aggfunc='nunique'),
    transactions_outside_cze=pd.NamedAgg(column='payment_outside_cze', aggfunc='mean'),
    transactions_inside_cze=pd.NamedAgg(column='payment_inside_cze', aggfunc='mean'),
    transactions_divisible_by_100=pd.NamedAgg(column='payment_divisible_by_100', aggfunc='sum'),
    total_transactions=pd.NamedAgg(column='amt', aggfunc='size'),
    average_distance=pd.NamedAgg(column='dist', aggfunc='mean'),
    average_hour_of_transaction=pd.NamedAgg(column='hour', aggfunc='mean')
)

grouped_df_by_stations = grouped_df_by_stations.rename(columns={
    'total_amount': 'celková částka',
    'average_transaction': 'průměrná transakce',
    'number_of_clients': 'počet unikátních klientů',
    'transactions_outside_cze': 'mimo ČR',
    'transactions_inside_cze': 'v ČR',
    'transactions_divisible_by_100': 'transakce delitelné 100',
    'total_transactions': 'celkový počet transakcí',
    'average_distance': 'průměrná vzdálenost cesty k ní',
    'average_hour_of_transaction': 'průměrná hodina transakce'
})

grouped_corr_matrix = grouped_df_by_stations.corr()

plt.figure(figsize=(16, 9))
sns.heatmap(grouped_corr_matrix, annot=True, cmap='coolwarm', fmt=".2f", annot_kws={"size": 10}, xticklabels=True, yticklabels=True)
plt.title('Korelační matice dat z benzínek', fontsize=20)
plt.xticks(fontsize=14, rotation=45, ha='right')
plt.yticks(fontsize=14)  
plt.show()

Zdá se, že samoobslužné beznínky jsou opravdu velice populární mezi klienty této banky.

Code

grouped_df_by_stations = df.groupby('posid').agg(
    number_of_payments=pd.NamedAgg(column='amt', aggfunc='size'),
    number_of_distinct_clients=pd.NamedAgg(column='clid', aggfunc='nunique'),
    average_amount=pd.NamedAgg(column='amt', aggfunc='mean')
)

plt.figure(figsize=(16, 8))

sns.histplot(grouped_df_by_stations[grouped_df_by_stations['average_amount'] < 2000]["average_amount"], kde=True, bins=30)
plt.title('Distribuce stanic dle průměrné výše platby', fontsize=20)
plt.xlabel('Průměrná částka', fontsize=16)
plt.ylabel('Počet stanic', fontsize=16)

plt.tight_layout()
plt.show()

Většina beznínek byla nevštívena pouze jedním jediným klientem či dvěmi.

Code

unique_clients_per_station = df.groupby('posid')['clid'].nunique()

bins = [0, 2, 5, 10, 15, 30, 70, 100, 200, 500, 1000, float('inf')]
labels = [f'{bins[i]+1} - {bins[i+1]}' for i in range(len(bins)-1)]

binned_data = pd.cut(unique_clients_per_station, bins=bins, labels=labels)

binned_count = binned_data.value_counts().sort_index()

binned_count_df = binned_count.reset_index()
binned_count_df.columns = ['Počet unikátních klientů', 'Počet stanic']
binned_count_df['Procenta'] = (binned_count_df['Počet stanic'] / binned_count_df['Počet stanic'].sum()) * 100

binned_count_df

	Počet unikátních klientů	Počet stanic	Procenta
0	1 - 2	16603	69.60
1	3 - 5	2524	10.58
2	6 - 10	1113	4.67
3	11 - 15	497	2.08
4	16 - 30	560	2.35
5	31 - 70	540	2.26
6	71 - 100	255	1.07
7	101 - 200	660	2.77
8	201 - 500	801	3.36
9	501 - 1000	262	1.10
10	1001 - inf	41	0.17

Podívejme se jaké benzínky jsou nejoblíbenější.

Code

df['name'] = df['name'].str.lower().str.strip()
text = " ".join(name for name in df.name.dropna() if ("cerpaci" not in name and "cs" not in name and "stanice" not in name))

wordcloud = WordCloud(collocations=False, width=800, height=400, background_color='white').generate(text)
plt.figure(figsize=(14, 7))
plt.imshow(wordcloud, interpolation='bilinear')
plt.axis("off")
plt.show()

Code

import numpy as np
total_amt_per_station = df.groupby(['pos_gps_lat', 'pos_gps_lon'])['amt'].sum().reset_index()

plt.figure(figsize=(14, 8))
scatter = plt.scatter(total_amt_per_station['pos_gps_lon'], total_amt_per_station['pos_gps_lat'], 
                      c=np.log(total_amt_per_station['amt']), cmap='hot', s=40, alpha=0.7)
plt.colorbar(scatter, label='Celkově utraceno na pumpách')
plt.title('Stanice dle lokace a celkové částky utracené na nich', fontsize=16)
plt.xlabel('Lon')
plt.ylabel('Lat')
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.show()

Časové série

Časové série umožnují získat vhled pro dlouhodobé trendy, cyklické vzory či nám pomoci identifikovat outliers a pak i případně dělat předpovědi pro dalši období (to ale udělám až v explanatorní části).

Zde jsou časové série jak vypadala průměrná transakce během roku a kolik jich bylo uděláno.

Code

from scipy.signal import find_peaks
df['date'] = pd.to_datetime(df['day'])

daily_totals = df.groupby(df['date'].dt.date)['amt'].sum()

daily_average = df.groupby(df['date'].dt.date)['amt'].mean()

daily_transactions = df.groupby(df['date'].dt.date).size()

def annotate_peaks(series, ax, height=50, prominence=500):
    peaks, _ = find_peaks(series.values, height=height, prominence=prominence)
    peak_dates = series.iloc[peaks].index
    for date in peak_dates:
        ax.annotate(date.strftime('%m-%d'), (date, series[date]), textcoords="offset points", xytext=(0,10), ha='center', fontsize=12)

def annotate_valleys(series, ax, height=50, prominence=500):
    inverted = series.max() - series
    valleys, _ = find_peaks(inverted.values, height=height, prominence=prominence)
    valley_dates = series.iloc[valleys].index
    for date in valley_dates:
        ax.annotate(date.strftime('%m-%d'), (date, series[date]), textcoords="offset points", xytext=(0,-15), ha='center', fontsize=12, color='blue')

plt.figure(figsize=(18, 15))

plt.subplot(2, 1, 1)
daily_average.plot(kind='line', color='green')
plt.title('Průměrná výše transakce v průběhu času', fontsize=15)
plt.ylabel('Průměrná částka', fontsize=13)
plt.xlabel('Datum', fontsize=13)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  
annotate_peaks(daily_average, plt.gca(), height=10, prominence=50)
annotate_valleys(daily_average, plt.gca(), height=10, prominence=50)


plt.subplot(2, 1, 2)
daily_transactions.plot(kind='line', color='red')
plt.title('Počet transakcí v průběhu času', fontsize=15)
plt.ylabel('Počet transakcí', fontsize=13)
plt.xlabel('Datum', fontsize=13)
plt.xticks(fontsize=14) 
plt.yticks(fontsize=14)  
annotate_peaks(daily_transactions, plt.gca(), height=500, prominence=500)
annotate_valleys(daily_transactions, plt.gca(), height=2500, prominence=700)


plt.tight_layout()
plt.show()

V časové sérii počtu transakcí (červený graf) je vidět jasná periodicita peaků, kterou lze odečíst z datumů u nich - jde většinou o pátky, kdy lidé jedou někam na víkend (jen opakování již vyobrazeného grafu).

Interpretovat peaky v prvním (zeleném grafu) zobrazujícím průměrnou částku je složitější. Peaky odpovídají jinému ,již předtím viděnému grafu, jenž zobrazoval počet a výše plateb dle dnu v měsíci. Nejspíše lidé utratí nejvíce po výplatě, která touto dobou chodí. Graf mi ale také připadá mnohem více zašuměný a tak dělat nějaké jasné závěry člověk úplně asi nemůže (minimálně já si na to netroufám).

Poklesy na úplných koncích jsou způsobeny Vánoci a oslavou nového roku. Peak 13 dubna 2017 přisuzuji tomu, že se jednalo o zelený čtvrtek a lidé tankovali před Velikonoci.

Explanatorní analýza

Clustering

Podívejme se jestli klienti vytvářejí nějaké clustery a jak vypadá vůbec varieta jejich chování. Z toho můžeme případně vyčíst jaké jsou typy klientů a jak si je jejich chování podobné.

Pro každého klienta byl vytvořen vektor těchto featur:

celková útrata, průměrná útrata, počet navštívených unikátních stanic a také celkově, počet navštívených unikátních zemí a opět i celkově, průmerná uražená vzdálenost a celkově uražená vzdálenost, modus hodiny kdy tankoval

Jelikož je dimenze vektoru 9, nelze ho nromálně vyzualizovat a je třeba provést dimensionální redukci - pro to jsem zvolil algoritmus PaCMAP.

Code

import pacmap
from sklearn.preprocessing import StandardScaler
import umap

df['dist'].fillna(0, inplace=True)
df.dropna(subset=['hour'], inplace=True)
grouped_df = df.groupby('clid').agg(
    total_spent=pd.NamedAgg(column='amt', aggfunc='sum'),
    avg_spent=pd.NamedAgg(column='amt', aggfunc='mean'),
    unique_stations=pd.NamedAgg(column='posid', aggfunc='nunique'),
    unique_countries=pd.NamedAgg(column='country', aggfunc='nunique'),
    total_stations_visited=pd.NamedAgg(column='posid', aggfunc='count'),
    total_countries_visited=pd.NamedAgg(column='country', aggfunc='count'),
    avg_distance=pd.NamedAgg(column='dist', aggfunc='mean'),
    total_distance=pd.NamedAgg(column='dist', aggfunc='sum'),
    median_hour=pd.NamedAgg(column='hour', aggfunc=lambda x: x.mode()[0] if not x.mode().empty else np.nan), #named median but is now mode
)

recalc = True
do_umap = False
encode_hours = False

# One-hot encode transaction hours and sum them up
if encode_hours:
    hour_dummies = pd.get_dummies(df['hour'], prefix='hour')
    hour_dummies['clid'] = df['clid']
    hour_sum = hour_dummies.groupby('clid').sum()

    # Join the calculated metrics and one-hot encoded hours
    final_df = pd.merge(grouped_df, hour_sum, on='clid')

else:
    final_df = grouped_df

# Standardize the features
scaler = StandardScaler()
scaled_df = scaler.fit_transform(final_df)

# Perform UMAP dimensionality reduction

if recalc:
    if do_umap:
        umap = umap.UMAP(n_neighbors=20, min_dist=0.1, n_components=2)
        umap_result = umap.fit_transform(scaled_df)
        np.save('umap_result.npy', umap_result)
    else:
        pacmap_result = pacmap.PaCMAP(n_components=2, n_neighbors=None, MN_ratio=0.5, FP_ratio=2.0).fit_transform(scaled_df)
        np.save('pacmap_result.npy', pacmap_result)
else:
    if do_umap:
        umap_result = np.load('umap_result.npy')
    else:
        pacmap_result = np.load('pacmap_result.npy')

if not do_umap:
    umap_result = pacmap_result

Předpokladem pro užití UMAPU a i PaCMAPu j, že body jsou na n-dimensionální Reimannovské varietě. Je otázka jak je tento předpoklad splněn (nikdy člověk tohohle ani nemá záruku). Předpokládám, že tady půjde spíše o pseudo-varietu, což ovšem moc nevadí - navíc UMAP i PaCMAP fungují dobře, i když je předpoklad porušen.

Code

umap_df = pd.DataFrame(umap_result, columns=['UMAP1', 'UMAP2'])
umap_df['clid'] = final_df.index
final_umap_df = pd.merge(umap_df, final_df.reset_index(), on='clid')
bins_avg_spent = [0, 50, 100, 200, 400, 700, 1100, 1300, np.inf]
amount = 30
bins_total_spent = [int(1.5**i) for i in range(10, amount)] + [np.inf]
bins_total_spent_labels = [f'{bins_total_spent[i]+1} - {bins_total_spent[i+1]}' for i in range(len(bins_total_spent)-1)]
final_umap_df['avg_spent_binned'] = pd.cut(final_umap_df['avg_spent'], bins_avg_spent, labels=['1-50', '50-100', '100-200','200-400', '400-700', '700-1100', '1100-1300', '1300+'])
final_umap_df['total_spent_binned'] = pd.cut(final_umap_df['total_spent'], bins_total_spent, labels=bins_total_spent_labels)
# print(final_umap_df['avg_spent_binned'])

fig, axs = plt.subplots(1, 2, figsize=(24, 10))

sns.scatterplot(x='UMAP1', y='UMAP2', hue='avg_spent_binned', data=final_umap_df, palette='viridis', legend="full", s = 10, linewidth=0.1, alpha=0.7, ax=axs[0])
axs[0].set_title('Colored by Average Amount Spent')
axs[0].legend(title='Average Spent (CZK)')

sns.scatterplot(x='UMAP1', y='UMAP2', hue='total_spent_binned', data=final_umap_df, palette='viridis', legend="full", s = 10, linewidth=0.1, alpha=0.7, ax=axs[1])
axs[1].set_title('Colored by Total Amount Spent')
axs[1].legend(title='Total Spent (CZK)')
plt.show()

Code

amount = 20
amount_2 = 7
bins_unique_stations = [i*3 for i in range(amount)] + [np.inf]  # Adjust bins based on your labels
bins_unique_stations_labels = [f"{bins_unique_stations[i]} - {bins_unique_stations[i+1]}" for i in range(len(bins_unique_stations)-1)]
bins_unique_countries = [i for i in range(amount_2)] + [np.inf]  # Adjust bins based on your labels
bins_unique_countries_labels = [f"{bins_unique_countries[i]} - {bins_unique_countries[i+1]}" for i in range(len(bins_unique_countries)-1)]
final_umap_df['unique_stations_binned'] = pd.cut(final_umap_df['unique_stations'], bins_unique_stations, labels=bins_unique_stations_labels)
final_umap_df['unique_countries_binned'] = pd.cut(final_umap_df['unique_countries'], bins_unique_countries, labels=bins_unique_countries_labels)

fig, axs = plt.subplots(1, 2, figsize=(24, 10))
sns.scatterplot(x='UMAP1', y='UMAP2', hue='unique_stations_binned', data=final_umap_df, palette='viridis', legend="full", s = 10, linewidth=0.1, alpha=0.7, ax=axs[0])
axs[0].set_title('Colored by Unique Stations Count')
axs[0].legend(title='Unique Stations Count')

sns.scatterplot(x='UMAP1', y='UMAP2', hue='unique_countries_binned', data=final_umap_df, palette='viridis', legend="full", s = 10, linewidth=0.1, alpha=0.7, ax=axs[1])
axs[1].set_title('Colored by Unique Countries Count')
axs[1].legend(title='Unique Countries Count')

plt.show()

Vidíme, že “rozvětvení” variety je kvůli počtu navštívených zemí.

Code

amount = 20
bins_avg_distance = [int(1.3**i) for i in range(4, amount+4)] + [np.inf]
bins_avg_distance_labels = [f"{bins_avg_distance[i]} - {bins_avg_distance[i+1]}" for i in range(len(bins_avg_distance)-1)]
bins_total_distance = [int(1.5**i) for i in range(2, amount+3)] + [np.inf]
bins_total_distance_labels = [f"{bins_total_distance[i]} - {bins_total_distance[i+1]}" for i in range(len(bins_total_distance)-1)]

final_umap_df['avg_distance_binned'] = pd.cut(final_umap_df['avg_distance'], bins_avg_distance, labels=bins_avg_distance_labels)
final_umap_df['total_distance_binned'] = pd.cut(final_umap_df['total_distance'], bins_total_distance, labels=bins_total_distance_labels)

fig, axs = plt.subplots(1, 2, figsize=(24, 10))

sns.scatterplot(x='UMAP1', y='UMAP2', hue='avg_distance_binned', data=final_umap_df, palette='viridis', legend="full", s = 10, linewidth=0.1, alpha=0.7, ax=axs[0])
axs[0].set_title('Colored by Average Distance')
axs[0].legend(title='Average Distance')

sns.scatterplot(x='UMAP1', y='UMAP2', hue='total_distance_binned', data=final_umap_df, palette='viridis', legend="full", s = 10, linewidth=0.1, alpha=0.7, ax=axs[1])
axs[1].set_title('Colored by Total Distance')
axs[1].legend(title='Total Distance')

plt.show()

“Paprsky” na levé straně jsou způsobeny hodinami kdy převážně klient byl na pumpě (jen zde nezobrazuji pro to specificky barvený plot, protože není zvlášť zajímavý).

Je vidět z plotu, jak lidé co utratili nejvíce peněz jsou zároveň ti co převážně hodně cestují (navštívili hodně benzínek, také byli v zahraničí a urážejí větší vzdálenosti). Na toto již poukazovala korelační matice. Lze je také považovat za jakýsi “středobod” chování jak ukazuje varieta - to dává i smysl, jelikož chování ostaních skupin je podmnožina jejich chování, pač dělají vlastně všechno co jde (nakupují v cizích zemích, tankují na hodně benzínách, atd.)

Nevypadá to ale, že by existovali klienti s nějakým vyloženě nespecifickým chováním a vyráběli by tak nějaký zvlášť vyhraněný cluster. Je možné, že v realitě takový je, ale tyto data neobsahují dostatečně podrobné informace o chování klientů (další údaje k transakcím), aby se to dalo zjistit.

Hledat clustery - případně jejich počet jde pomocí počítaní takzvané gap statistiky.

Code

from sklearn.cluster import KMeans
from gap_statistic import OptimalK
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaled_data = scaler.fit_transform(grouped_df)

optimalK = OptimalK()
n_clusters = optimalK(scaled_data, cluster_array=np.arange(1, 100))
gap_df = optimalK.gap_df

plt.figure(figsize=(10, 6))
plt.plot(optimalK.gap_df.n_clusters, optimalK.gap_df.gap_value, linewidth=3)
plt.scatter(optimalK.gap_df[optimalK.gap_df.n_clusters == n_clusters].n_clusters,
            optimalK.gap_df[optimalK.gap_df.n_clusters == n_clusters].gap_value, s=200, c='r')
plt.grid(True)
plt.xlabel('Cluster Count')
plt.ylabel('Gap Value')
plt.title('Gap Values by Cluster Count')
plt.show()

Gap statistika zde také poukazuje na to, že neexistuje žádný zvlášť “dobrý” počet clusterů.

Tohle samé by pak šlo dělat pro benzínky, ale neočekávám, že bych našel něco zajimáveho a bylo by to podobné tomuto vzhledem k datové kvalitě. Například mi chybí vysoce relevatní údaje pro clustering jako jestli je beznínka ve městě, na vesnici nebo u dálnice - jaké nabízela ceny za litr v dobu čepování a další podobné. Nyní mám vlastně pouze kolik tam lidi utratili a kdy. Tohle by šlo doscrapovat z internetu (dle polohy a názvu), ale to by bylo mimo rozsah této práce, tak to pouze zmiňuji.

Predikce plateb

GBDT

Zde vyrobím model, který bude předpovídat, kolik celkově utratí člověk daný měsíc, pokud víme chování klienta v daném měsíci (kolik navštívil pump, v jakých zemích atd.), jen nevíme kolik utrácel.

Nejprve je třeba data opět transformovat a vyrobit si další featury.

Code

from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error
monthly_payments = df.groupby(['clid', 'month']).agg(total_spent=pd.NamedAgg(column='amt', aggfunc='sum'),
                                                        unique_stations=pd.NamedAgg(column='posid', aggfunc='nunique'),
                                                        unique_countries=pd.NamedAgg(column='country', aggfunc='nunique'),
                                                        total_stations_visited=pd.NamedAgg(column='posid', aggfunc='count'),
                                                        total_countries_visited=pd.NamedAgg(column='country', aggfunc='count'),
                                                        avg_distance=pd.NamedAgg(column='dist', aggfunc='mean'),
                                                        total_distance=pd.NamedAgg(column='dist', aggfunc='sum'),
                                                        mod_hour=pd.NamedAgg(column='hour', aggfunc=lambda x: x.mode()[0] if not x.mode().empty else np.nan),
                                                     ).reset_index()

pivot_df = monthly_payments.pivot_table(index='clid', columns='month', values='total_spent')

pivot_df.columns = ['total_spent_' + str(col) for col in pivot_df.columns]
pivot_df = pivot_df.reset_index()

monthly_payments = pd.merge(monthly_payments, pivot_df, on='clid', how='left')
total_spent_columns = [col for col in monthly_payments.columns if 'total_spent_' in col]

#avoid data leakage
for row in monthly_payments.itertuples():
    for col in total_spent_columns:
        if "total_spent_"+str(row.month) == col:
            monthly_payments.at[row.Index, col] = np.nan

monthly_payments['avg_spent_other_months'] = monthly_payments[total_spent_columns].apply(lambda row: row.mean(), axis=1)
monthly_payments['total_spent_other_months'] = monthly_payments[total_spent_columns].apply(lambda row: row.sum(), axis=1)
monthly_payments = pd.get_dummies(monthly_payments, columns=['mod_hour', "month"])

Takto vypadá dataframe, jenž bude užit pro trenování s targetem: kolik bylo utraceno daný měsíc. Nans byly nahrazeny hodnotou -1. Měsíce a hodiny byly onehot enkódovány. A pak byla data normalizována.

Code

features = monthly_payments.columns.tolist()
features.remove('total_spent')
features.remove('clid')

# monthly_payments = monthly_payments[monthly_payments['total_spent'] != -1]
monthly_payments.fillna(-1, inplace=True)
monthly_payments_scaled = pd.DataFrame(scaler.fit_transform(monthly_payments), columns=monthly_payments.columns)

X = monthly_payments_scaled[features]
y = monthly_payments['total_spent']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
X_train

	unique_stations	unique_countries	total_stations_visited	total_countries_visited	avg_distance	total_distance	total_spent_1	total_spent_2	total_spent_3	total_spent_4	...	month_3	month_4	month_5	month_6	month_7	month_8	month_9	month_10	month_11	month_12
165402	-0.59	-0.18	0.76	0.76	-0.54	-0.37	0.37	-0.08	1.42	0.12	...	-0.28	-0.29	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	-0.31	3.23
597878	-0.59	-0.18	-0.60	-0.60	0.62	-0.12	-0.58	-0.61	-0.65	-0.65	...	-0.28	-0.29	-0.30	-0.30	3.19	-0.32	-0.31	-0.31	-0.31	-0.31
74417	-0.59	-0.18	0.76	0.76	-0.57	-0.40	1.37	0.41	0.02	2.22	...	-0.28	-0.29	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	3.25	-0.31
38563	1.31	-0.18	0.76	0.76	0.25	0.59	0.95	-0.12	-0.65	0.34	...	-0.28	-0.29	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	-0.31	3.23
661924	-0.59	-0.18	-0.60	-0.60	-0.07	-0.29	-0.69	-0.61	-0.65	-0.68	...	-0.28	-0.29	-0.30	-0.30	-0.31	-0.32	3.21	-0.31	-0.31	-0.31
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
644167	-0.59	-0.18	-0.60	-0.60	-0.45	-0.38	-0.69	-0.61	-0.65	-0.68	...	-0.28	-0.29	-0.30	3.28	-0.31	-0.32	-0.31	-0.31	-0.31	-0.31
259178	0.04	-0.18	0.08	0.08	0.20	0.15	-0.63	0.31	1.01	-0.68	...	-0.28	3.46	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	-0.31	-0.31
365838	-0.59	-0.18	-0.26	-0.26	-0.55	-0.39	-0.18	-0.61	-0.65	-0.05	...	3.56	-0.29	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	-0.31	-0.31
131932	-0.59	-0.18	-0.26	-0.26	-0.51	-0.38	-0.69	0.67	-0.65	0.62	...	-0.28	-0.29	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	-0.31	-0.31
121958	-0.59	-0.18	-0.60	-0.60	-0.57	-0.41	-0.69	-0.17	1.69	-0.68	...	-0.28	-0.29	-0.30	-0.30	-0.31	-0.32	-0.31	-0.31	3.25	-0.31

534500 rows × 56 columns

Jako model jsem zde zvolil Gradient Boosted Regression Decision Tree, protože se ukazuje jako nejlepší model pro tabulární data a v mém testování i přinesl ty nejelpší výsledky a porazil náhodné lesy či SVM.

Pro vyhodnocení modelu jsem zvolil Root Mean Absolute Error vzhledem k typu problému a dobré intepretovatelnosti. Poměr train a test setu jsem zvolil 80:20.

Každopádně bohužel model nefunguje příliš přesně a chyba je poměrně vysoká. Vidím za tím primárně dva důvody: Nedostatek featur - ono si můžu přidělat více featur (statistiky pro předchozí měsíce ale to je pak 90+ featur), ale můj počítač to trénoval už pak staršně dlouho (už tohle trvalo tako 2O minut a přece jenom tady jde o demonstraci úlohy…). Druhým důvodem je, že data jsou velice šuměná, což je primárně způsobeno náhodou.

Code

import pickle
import os
import lzma
if os.path.exists('gbm_client_pred_model.pkl'):
    with lzma.open('gbm_client_pred_model.pkl', 'rb') as f:
        model = pickle.load(f)
else:
    model = GradientBoostingRegressor(n_estimators=100, max_depth=6, min_samples_split=2, learning_rate=0.1)
    model.fit(X_train, y_train)
    with lzma.open('gbm_client_pred_model.pkl', 'wb') as f:
        pickle.dump(model, f)

predictions = model.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, predictions))
print(f'RMSE: {rmse}')


df_pred = pd.DataFrame({'Predictions': predictions, 'Actual': y_test})

random_subset = df_pred.sample(n=10)

random_subset

RMSE: 868.306379195834

	Predictions	Actual
41270	1507.56	1564.25
89959	864.66	584.10
421462	1643.53	2000.20
188481	866.35	366.20
62584	649.09	378.39
220181	3631.11	4464.02
348865	2925.14	2831.45
588740	1413.34	2000
160310	1916.03	2900.10
356073	283.53	300

Překvapivě jedna z nejdůležitějších featur pro predikci je to kolikrát klient pumpoval v cizině.

Code

total_spent_weights = 0
mod_hour_weights = 0
month_weights = 0
for feature, importance in zip(features, model.feature_importances_):
    if "total_spent" in feature:
        total_spent_weights += importance
    elif "mod_hour" in feature:
        mod_hour_weights += importance
    elif "month" in feature:
        month_weights += importance
    else:
        print(f'{feature}: {importance:.4f}')
print(f'Total spent weights for each month: {total_spent_weights:.4f}')
print(f'Mod hour weights for each month: {mod_hour_weights:.4f}')
print(f'Month weights for each month: {month_weights:.4f}')

unique_stations: 0.0178
unique_countries: 0.0243
total_stations_visited: 0.0846
total_countries_visited: 0.3407
avg_distance: 0.0027
total_distance: 0.0078
Total spent weights for each month: 0.0628
Mod hour weights for each month: 0.0056
Month weights for each month: 0.4536

Knn

O něco praktičtější je zkusit predikovat kolik daný klient utratí v daný měsíc, když víme jeho chování v ostatních měsícíh. Ovšem kdybych chtěl dělat model, který jen a pouze dostane na vstupu klienta a předpovídal bych kolik utratí v daný měsíc, tak výsledek bude ještě horší než u předchozího problému inherentně (protože uvidím ještě méně dat). Přesto ale model může být užitečný - zvolím přístup pomocí Knn a najdu lidi, kteří vykazují blízké chování danému klientovi a z jejich utráty pro daný měsíc předpovím to jak utrácel náš klient.

Code

from sklearn.neighbors import KNeighborsRegressor

monthly_df = df.groupby(['clid', 'month']).agg(total_spent=pd.NamedAgg(column='amt', aggfunc='sum'),
                                                        unique_stations=pd.NamedAgg(column='posid', aggfunc='nunique'),
                                                        unique_countries=pd.NamedAgg(column='country', aggfunc='nunique'),
                                                        total_stations_visited=pd.NamedAgg(column='posid', aggfunc='count'),
                                                        total_countries_visited=pd.NamedAgg(column='country', aggfunc='count'),
                                                        avg_distance=pd.NamedAgg(column='dist', aggfunc='mean'),
                                                        total_distance=pd.NamedAgg(column='dist', aggfunc='sum'),
                                                        mod_hour=pd.NamedAgg(column='hour', aggfunc=lambda x: x.mode()[0] if not x.mode().empty else np.nan),
                                                     ).reset_index()
# Pivot the DataFrame
pivot_df = monthly_df.pivot(index='clid', columns='month')
pivot_df = pivot_df.reset_index()
pivot_df.columns = ['_'.join(str(col) for col in col_tuple).strip() for col_tuple in pivot_df.columns.values]

Trénovací data jsou podobná těm předchozím, ale obsahují veškeré info z předchozích měsíců (tady si to můžu dovolit, protože to bude rychlé i tak) a vektor pro jednoho člověka je tedy 88 velký. Zde predikuji data pro červen (jsem si náhodně vybral tento měsíc).

Code

prediction_month = "6" # Change this to the month you want to predict

scaler = StandardScaler()
pivot_df.fillna(-1, inplace=True)
# pivot_df = pivot_df[pivot_df['total_spent_'+prediction_month] != -1]
pivot_df_scaled = pd.DataFrame(scaler.fit_transform(pivot_df), columns=pivot_df.columns)

cols_to_include = [col for col in pivot_df.columns if prediction_month not in col and "total_spent" in col]
# cols_to_include.remove('clid_')
X = pivot_df_scaled.loc[:, cols_to_include]
y = pivot_df['total_spent_'+prediction_month]  # Target month spending
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
X_train

	total_spent_1	total_spent_2	total_spent_3	total_spent_4	total_spent_5	total_spent_7	total_spent_8	total_spent_9	total_spent_10	total_spent_11	total_spent_12
22982	-0.55	-0.47	-0.18	-0.52	-0.53	-0.56	-0.55	-0.56	-0.56	-0.55	-0.57
90478	-0.53	-0.47	-0.50	-0.52	-0.53	-0.56	-0.54	-0.56	-0.56	-0.55	-0.57
49685	-0.55	-0.47	-0.50	-0.52	-0.53	-0.56	-0.41	-0.56	-0.56	-0.55	-0.57
114701	-0.55	-0.47	-0.50	-0.52	-0.53	-0.56	-0.55	-0.31	-0.25	-0.55	-0.57
47661	-0.55	0.70	-0.50	-0.11	-0.53	-0.56	-0.55	0.61	-0.17	-0.49	-0.57
...	...	...	...	...	...	...	...	...	...	...	...
110268	-0.55	-0.47	-0.50	-0.52	-0.19	-0.56	-0.55	-0.48	-0.56	-0.08	-0.57
119879	-0.49	-0.47	-0.50	-0.52	-0.53	-0.56	-0.55	-0.56	-0.43	-0.55	-0.57
103694	-0.55	-0.47	-0.50	-0.27	-0.53	-0.56	-0.55	-0.56	-0.56	-0.24	-0.57
860	-0.55	-0.47	-0.50	-0.52	0.06	-0.56	-0.55	-0.56	-0.56	-0.55	-0.57
15795	0.27	2.93	1.60	4.08	2.27	2.45	1.71	1.21	1.85	1.29	1.13

97431 rows × 11 columns

Hyperparametr pro počet sousedů jsem zvolil 10.

Code

knn = KNeighborsRegressor(n_neighbors=10, weights="uniform")
knn.fit(X_train, y_train)

KNeighborsRegressor(n_neighbors=10)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Code

y_pred = knn.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
print(f'Root Mean Squared Error: {np.sqrt(mse)}')

Root Mean Squared Error: 871.8034504118918

Překvapivě síla modelu je poměrně podobně kvalitná jako i GBDT. Důležité je ale poznamenat, že zde výsledky nejsou úplně srovnatelné, protože GBDT predikuje utrátu náhodně ze všech měsíců, zatímco tady jsem zvolil pouze červen, ale lze čekat, že výsledky by byly podobné i kdyby dělaly ty samé predikce.

Případný neduh tohoto modelu a i toho předchozího je ten, že rozumně nepředpoví Nan, tedy to když klient nic nenakoupil, ovšem reálně asi nejde o takový problém, protože ho budeme chtít využít pravě pro klienty, kde čekáme, že daný měsíc nějak něco utráceli.

Baseline

Dobré je si výsledky modelů srovnat s trivialním baseline modelem, který vždy předpoví průměrnou hodnotu.

Code

X_for_mean = pivot_df.loc[:, cols_to_include]
y_for_mean = pivot_df['total_spent_'+prediction_month]

X_for_mean_train, X_for_mean_test, y_for_mean_train, y_for_mean_test = train_test_split(X_for_mean, y_for_mean, test_size=0.2, random_state=42)

y_mean_pred = X_for_mean_test.mean(axis=1)
mse = mean_squared_error(y_for_mean_test, y_mean_pred)
print(f'Root Mean Squared Error: {np.sqrt(mse)}')

Root Mean Squared Error: 864.4262627749122

Ukazuje se, že baseline má dokonce lepší výsledky než Knn. Doplnění hodnoty náhodného jevu (utracená částka za měsíc) pomocí průměru je opravdu konzistetní způsob - ostatně, je to nejvíce konzistetní odhad dle MLE. Ovšem zde jevy nejsou úplně nezavávislé a tedy MLE předpoklad je porušen - důvodem je například i to, že v některé měsíce lidé tankují za více, protože je zima, jak již ukázali grafy v exploratorní analýze, i tak ale je baseline model stejně výkonný jak mnou udělané zbylé výše.

Závěr

Dataset je plný informací pro zjištění a dokázání nějakých základních vzorců chování lidí na beznínových pumpách - například to, že v zimě se utrací na beznínce v průměru více než v jiné roční období, že lidé tankují nejvíce v pátek nebo to, že většina plateb je do pětiset korun a také to kdy lidé převážně tankují a zakolik. Téže bylo nalezeno možné spojení mezi samoobslužnými beznínkami a hezkými částkami a problémy v platbách (duplicity či chybějící hodiny) nebo důvod růstu či poklesu držeb v dané dny kvůli svátkům.

Pokročilé modely se ukázaly ne jako příliš praktické, protože lidé autem jezdí poměrně stejně každý měsíc a navíc je v datech šum například, že zaplatí na pumpě manželka a já pak vidím pokles v tom kolik utratil klient daný měsíc, protože si tohle nemám jak propojit. Proto se také odhad pomocí průměru ukázal jako dobrý prediktor kolik člověk utratí/l daný měsíc.

Bohužel v datasetu chybí zajímavé a podstatné informace, které by mohli beznínovým pumpám napomoci k optimalizaci chování a tím zvýšení zisku. Například to jaká byla cena benzínu za litr v době tankování, jaká je jeho lokolita (dálnice, město…). Také chybí informace o tom, jaké služby beznínky nabízí (například samoobslužné pumpy, či jiné) a jaké položky byly nakoupeny v transakci. Něco z tohohle by šlo samozřejmě doplnit z internetu, ale to by bylo mimo rozsah této práce jak jsem již avizoval. Myslím, ale že by tato práce mohla sloužit jako dobrý základ pro další analýzy, které by se mohly zaměřit na další aspekty.