<jupyter_start><jupyter_text>DATA PREPARATION<jupyter_code>df = pd.read_csv('./Homework/Data/finalData.csv')
df.info()
print('Transactions were made from {} to {}'.format(df['OrderDate'].min(), df['OrderDate'].max()))
print("Result of distinct value:\n", df.nunique())
df['TotalValue'] = df['OrderQty']*(1 - df['UnitPriceDiscount'])*df['UnitPrice']
(df['TotalValue'] - df['LineTotal']).unique()<jupyter_output><empty_output><jupyter_text>Because the subtraction between each row of TotalValue and LineTotal is almost 0, we can assume that LineTotal contains final value of transactions, which will be used to calculate Monetary factor later.<jupyter_code>df.drop('TotalValue', axis=1, inplace=True)
df['OrderDate'] = pd.to_datetime(df['OrderDate'], dayfirst=True)
df['OrderDate'].describe()
df[df['CustomerID']==11000]<jupyter_output><empty_output><jupyter_text>For a random customer, we see that if he/she bought more than one product, the SalesOrderID and OrderDate will be duplicated for each ProductID. However, when counting Frequency, we still recognize that this customer made 3 unique purchases.<jupyter_code>datemark = dt.datetime(2021,12,26)
df['Recency'] = df['OrderDate']
rfm = df.groupby('CustomerID').agg({'Recency': lambda date: (datemark - date.max()).days,
                                    'T': lambda date: (datemark - date.min()).days,
                                    'SalesOrderID': lambda num: num.nunique(),
                                    'LineTotal': lambda price: price.sum()}).reset_index()
rfm = rfm.rename(columns={'SalesOrderID': 'Frequency', 'LineTotal': 'Monetary'})
rfm<jupyter_output><empty_output><jupyter_text>TRADITIONAL RFM MODEL   EDA<jupyter_code>rfm.describe()
print("Result of duplicated: ", rfm.duplicated().sum())
print("-------------")
print("Result of missing value: ", rfm.isnull().values.sum())
print("-------------")
print("Result of distinct value:\n", rfm.nunique())
print("-------------")

#boxplot
list_num =['Recency','Frequency','Monetary']
for col in list_num:
    plt.boxplot(rfm[col], vert=False, notch=True, flierprops=dict(markerfacecolor='orange', marker='s'))
    plt.title(col)
    plt.show()

for column in ['Recency', 'Frequency', 'Monetary']:
    print(f"Sample Variance: {np.var(rfm[column])}")
    
    print(f"SE:", np.std(rfm[column], ddof=1) / np.sqrt(np.size(rfm[column])))
    print(f"Mode: {mode(rfm[column])}")
    print(f"Median: {sta.median(rfm[column])}")

    print(f"Skewness: {skew(rfm[column], bias = False)}")
    print(f"Kurtosis: {kurtosis(rfm[column], bias = False)}")
    print("--------------")
for i in ['Recency','Frequency','Monetary']:
    q25, q75 = np.percentile(rfm[i], 25), np.percentile(rfm[i], 75)
    iqr = q75 - q25
    print('Percentiles: 25th=%.3f, 75th=%.3f, IQR=%.3f' % (q25, q75, iqr))
    cut_off = iqr * 1.5
    lower, upper = q25 - cut_off, q75 + cut_off
    # identify outliers
    outliers = [x for x in rfm[i] if x < lower or x > upper]
    print('Identified outliers: %d' % len(outliers))
    # non-outliers
    non_outliers = [x for x in rfm[i] if x >= lower and x <= upper]
    print('Non-outlier observations: %d' % len(non_outliers))<jupyter_output><empty_output><jupyter_text>The percentage of outliers in each column is in range between 2.4% and 4.4%. However, we decide to keep them first.<jupyter_code>quantiles = rfm[['Recency', 'Frequency', 'Monetary']].quantile(q=[0.2, 0.4, 0.6, 0.8]).to_dict()
quantiles
def r_score(x):
    if x <= quantiles['Recency'][.2]:
        return 5
    elif x <= quantiles['Recency'][.4]:
        return 4
    elif x <= quantiles['Recency'][.6]:
        return 3
    elif x <= quantiles['Recency'][.8]:
        return 2
    else:
        return 1

def fm_score(x, c):
    if x <= quantiles[c][.2]:
        return 1
    elif x <= quantiles[c][.4]:
        return 2
    elif x <= quantiles[c][.6]:
        return 3
    elif x <= quantiles[c][.8]:
        return 4
    else:
        return 5    

rfm['R'] = rfm['Recency'].apply(lambda x: r_score(x))
rfm['F'] = rfm['Frequency'].apply(lambda x: fm_score(x, 'Frequency'))
rfm['M'] = rfm['Monetary'].apply(lambda x: fm_score(x, 'Monetary'))
rfm['RFM_Segment'] = rfm['R'].map(str) + rfm['F'].map(str) + rfm['M'].map(str)
rfm['RFM_Score'] = rfm[['R', 'F', 'M']].sum(axis=1)
rfm.groupby('RFM_Segment').agg({'Recency': 'mean', 'Frequency': 'mean', 'Monetary': ['mean', 'count']}).round(1)
segt_map = {
    r'[1-2][1-2]': 'hibernating',
    r'[1-2][3-4]': 'at risk',
    r'[1-2]5': 'can\'t loose',
    r'3[1-2]': 'about to sleep',
    r'33': 'need attention',
    r'[3-4][4-5]': 'loyal customers',
    r'41': 'promising',
    r'51': 'new customers',
    r'[4-5][2-3]': 'potential loyalists',
    r'5[4-5]': 'champions'
}

rfm['Segment'] = rfm['R'].map(str) + rfm['F'].map(str)
rfm['Segment'] = rfm['Segment'].replace(segt_map, regex=True)
rfm_plot = rfm.groupby('Segment').agg(customers=('CustomerID', 'count'), min_rfm=('RFM_Segment', 'min'), max_rfm=('RFM_Segment', 'max')).sort_values(by='customers', ascending=False).reset_index()
rfm_plot
treemap = rfm.groupby('Segment').agg(customers=('CustomerID', 'count')).reset_index()
import squarify
from matplotlib.pyplot import figure
figure(figsize=(35, 15), dpi=300)
labels=treemap['Segment']
ax = squarify.plot(sizes=treemap['customers'], label=labels, alpha=1, text_kwargs={'color': 'black', 'size': 25},
              color=plt.cm.RdYlGn(np.linspace(0, 1, len(labels))), ec='black', norm_x=145, norm_y=90)
ax.axis('off')
bars = plt.barh(width=rfm_plot['customers'], y=list(rfm_plot['Segment'].values), color='silver')
for i, bar in enumerate(bars):
    if rfm['Segment'].value_counts().index[i] in ['loyal customers', 'champions']:
        bar.set_color('orange')
plt.title('Number of customers per segmentation')
plt.show()<jupyter_output><empty_output><jupyter_text>KMEANS CLUSTERING  DATA TRANSFORMATION<jupyter_code># Distribution before transformation
plt.figure(figsize=(9, 12))

plt.subplot(3, 1, 1)
plt.title('Distribution of Recency')
sns.distplot(rfm['Recency'], hist_kws=dict(edgecolor="k", linewidth=1.2))

plt.subplot(3, 1, 2)
plt.title('Distribution of Frecency')
sns.distplot(rfm['Frequency'], hist_kws=dict(edgecolor="k", linewidth=1.2))

plt.subplot(3, 1, 3)
plt.title('Distribution of Monetary')
sns.distplot(rfm['Monetary'], hist_kws=dict(edgecolor="k", linewidth=1.2))

plt.tight_layout()
# Transformation process to get the skewness close to 0
# Log transformation
rfm_log = rfm.copy()
rfm_log = np.log(rfm_log[['Recency','Frequency','Monetary']])
# Square root transformation
rfm_sqrt = rfm.copy()
rfm_sqrt = np.sqrt(rfm_sqrt[['Recency','Frequency','Monetary']])
# Cube root transformation
rfm_cbrt = rfm.copy()
rfm_cbrt = np.cbrt(rfm_cbrt[['Recency','Frequency','Monetary']])
def check_skew(df, column):
    skewness = skew(df[column])
    plt.title('Distribution of ' + column)
    sns.distplot(df[column], hist_kws=dict(edgecolor="k", linewidth=1.2))
    print("{}'s skew: {}".format(column, skewness))
    return
# Square root transformation
plt.figure(figsize=(9, 9))

plt.subplot(3, 1, 1)
check_skew(rfm_sqrt,'Recency')

plt.subplot(3, 1, 2)
check_skew(rfm_sqrt,'Frequency')

plt.subplot(3, 1, 3)
check_skew(rfm_sqrt,'Monetary')

plt.tight_layout()
# Cube root transformation
plt.figure(figsize=(9, 9))

plt.subplot(3, 1, 1)
check_skew(rfm_cbrt,'Recency')

plt.subplot(3, 1, 2)
check_skew(rfm_cbrt,'Frequency')

plt.subplot(3, 1, 3)
check_skew(rfm_cbrt,'Monetary')

plt.tight_layout()
# Log transformation
plt.figure(figsize=(9, 9))

plt.subplot(3, 1, 1)
check_skew(rfm_log,'Recency')

plt.subplot(3, 1, 2)
check_skew(rfm_log,'Frequency')

plt.subplot(3, 1, 3)
check_skew(rfm_log,'Monetary')

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>It is clear that log-transformation gives us the skewness closest to 0. As a result, log-transformation will be chosen as a method to transform our data.<jupyter_code># Normalize data
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(rfm_log)
rfm_scaled = scaler.transform(rfm_log)
rfm_scaled = pd.DataFrame(rfm_scaled, columns=rfm_log.columns)
rfm_scaled.describe()<jupyter_output><empty_output><jupyter_text>FIND THE OPTIMAL K<jupyter_code>from scipy.spatial.distance import cdist
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
inertias = [] 
for k in range(1, 12):
    km = KMeans(n_clusters=k, random_state=1)
    km.fit(rfm_scaled)
    inertias.append(km.inertia_) 
plt.figure(figsize=(9, 9))

plt.subplot(2, 1, 2)
plt.plot(range(1, 12), inertias, 'bx-') 
plt.xlabel('Values of K') 
plt.ylabel('Inertia') 
plt.title('The Elbow Method using Inertia') 
plt.tight_layout()
from yellowbrick.cluster import KElbowVisualizer
km = KMeans(random_state=42)
visualizer = KElbowVisualizer(km, k=(2, 12))
visualizer.fit(rfm_scaled)
visualizer.show()
from yellowbrick.cluster import SilhouetteVisualizer
fig, ax = plt.subplots(3, 2, figsize=(15, 8))
for i in range(6, 10):
    km = KMeans(n_clusters=i, init='k-means++', random_state=42)
    q, mod = divmod(i, 2)
    visualizer = SilhouetteVisualizer(km, colors='yellowbrick', ax=ax[q-1][mod])
    visualizer.fit(rfm_scaled)
sh_score = []
for i in range(2, 12):
    km = KMeans(n_clusters = i, random_state=1)
    km.fit_predict(rfm_scaled)
    score = silhouette_score(rfm_scaled, km.labels_, metric='euclidean')
    sh_score.append(score)

plt.plot(sh_score)
plt.xticks(np.arange(10), range(2, 12))<jupyter_output><empty_output><jupyter_text>VISUALIZE SEGMENTATION<jupyter_code>import plotly.express as px

def kmeans(df, clusters_num):
    km = KMeans(n_clusters = clusters_num, random_state=1)
    y_clusters = km.fit_predict(df)
    print(f"Silhouette Score for {clusters_num}:", "%1.2f" % silhouette_score(rfm_scaled, clusters))
    # print(kmean.labels_)
    # print(kmean.cluster_centers_)
    fig = px.scatter_3d(df, x="Recency", y="Frequency", z="Monetary", color = clusters + 1, size = clusters + 1)
    # color, size increment to 1 to have the calibration starting from 1, instead of zero
    fig.show()
plane = dict(xaxis = dict(title  = 'Recency -->'), yaxis = dict(title  = 'Frequency --->'),zaxis = dict(title  = 'Monetary-->'))
import plotly.graph_objects as go
def kmeans(df, clusters_num):
    km = KMeans(n_clusters = clusters_num, random_state=1)
    labels = km.labels_
    trace = go.Scatter3d(x=x[:, 0], y=x[:, 1], z=x[:, 2], mode='markers',marker=dict(color=labels, size=10, line=dict(color='black', width=10)))
    layout = go.Layout(margin=dict(l=0, r=0), scene=plane, height=800, width=800)
    data = [trace]
    fig = go.Figure(data = data, layout = layout)
    fig.show()
km = KMeans(n_clusters = 4, random_state=1)
km.fit_predict(rfm_scaled)
labels = km.labels_
trace = go.Scatter3d(x=x[:, 0], y=x[:, 1], z=x[:, 2], mode='markers',marker=dict(color=labels, size=10, line=dict(color='black', width=10)))
layout = go.Layout(margin=dict(l=0, r=0), scene=plane, height=800, width=800)
data = [trace]
fig = go.Figure(data = data, layout = layout)
fig.show()
kmeans(rfm_scaled, 6)<jupyter_output><empty_output><jupyter_text>CUSTOMER LIFETIME VALUE - The CLTV account has many calculation metrics. Basically, the metrics to focus on: Average Order Value, Purchase Frequency, Profit Margin, Churn Rate, Customer Value, Customer Life Time Value.- Formula:   $CLTV = (Customer Value  /  Churn Rate) * Profit Margin$      $Customer Value = Average Order Value * Purchase Frequency$      $Average Order Value = Total Revenue * Total Number of Orders$      $Purchase Frequency = Total Number of Orders  /  Total Number of Customers$      $Profit Margin = Total Price  /  Profit Rate$      $Churn Rate = 1 - Repeat Rate$      $Repeat Rate = Number of Customers Who have Shopped more than Once / Total Number of Customers$  COHORT ANALYSIS<jupyter_code>df['month'] = df['OrderDate'].dt.month
df['year'] = df['OrderDate'].dt.year
print(df)
df['cohort'] = df.apply(lambda row: row['year'] * 100 + row['month'], axis=1)
cohorts = df.groupby('CustomerID')['cohort'].min().reset_index()
print(cohorts)
cohorts.columns = ['CustomerID', 'first_cohort']
df = df.merge(cohorts, on='CustomerID', how='left')

headers = df['cohort'].value_counts().reset_index()
headers.columns = ['Cohorts', 'Count']
headers = headers.sort_values(['Cohorts'])['Cohorts'].to_list()
df['cohort_distance'] = df.apply(lambda row: (headers.index(row['cohort']) - headers.index(row['first_cohort'])) if (row['first_cohort'] != 0 and row['cohort'] != 0) else np.nan, axis=1)
cohort_pivot = pd.pivot_table(df, index='first_cohort', columns='cohort_distance', values='CustomerID', aggfunc=pd.Series.nunique)
cohort_pivot = cohort_pivot.div(cohort_pivot[0], axis=0)
cohort_pivot
fig, ax = plt.subplots(figsize=(30, 30))
y_labels = [str(int(header) % 100) + '-' + str(int(header) / 100) for header in headers]
x_labels = range(0, len(y_labels))
plt.yticks(ticks=headers, labels=y_labels)
plt.xticks(x_labels, x_labels)
ax.set(xlabel='Months After First Purchase', ylabel='First Purchase Cohort')
plt.title("Cohort Analysis", fontsize=20)
sns.heatmap(cohort_pivot, annot=True, fmt='.0%', mask=cohort_pivot.isnull(), ax=ax, square=True, linewidths=1, cmap=sns.cubehelix_palette(rot=-.4))

plt.show()<jupyter_output><empty_output><jupyter_text>PARETO/NBD COMBINED WITH GAMMA - GAMMA MODEL<jupyter_code>rfm[['Frequency', 'Monetary']].corr()<jupyter_output><empty_output><jupyter_text>The Gamma-Gamma model is based on the assumption that the number of transactions does not depend on their monetary value.The frequency and monetary value are not correlated if the output is close to zero.<jupyter_code>rfm = rfm[rfm['Monetary'] > 0]
from lifetimes import ParetoNBDFitter
pareto_nbd = ParetoNBDFitter(penalizer_coef = 0.0)
pareto_nbd.fit(rfm["Frequency"], rfm["Recency"], rfm["T"])
from lifetimes.plotting import plot_frequency_recency_matrix
from lifetimes.plotting import plot_probability_alive_matrix
from lifetimes.plotting import plot_period_transactions

plt.figure(figsize=(8,6))
plot_frequency_recency_matrix(pareto_nbd)
plot_probability_alive_matrix(pareto_nbd)
plot_period_transactions(pareto_nbd)
rfm["p_not_alive"] = 1-pareto_nbd.conditional_probability_alive(rfm["Frequency"], rfm["Recency"], rfm["T"])
rfm["p_alive"] = pareto_nbd.conditional_probability_alive(rfm["Frequency"], rfm["Recency"], rfm["T"])
rfm["predicted_purchases"] = pareto_nbd.conditional_expected_number_of_purchases_up_to_time(30, rfm["Frequency"], rfm["Recency"], rfm["T"])
rfm.sort_values(by = "predicted_purchases")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1704685996976776, 'fraction_numerical': 0.02047547066098667, 'mean_word_length': 5.491971454058876, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '237784', 'license': 'CC0-1.0', 'n_tokens_mistral': 5683, 'n_tokens_neox': 5368, 'n_words': 1123, 'path': 'Homework/CustomerSegmentation_CLV.ipynb', 'repo_name': 'tarynduong/projects-written-in-python', 'size': '14554'}