Customer Segmentation and Marketing with Markov Chains in Python: A Data Science Tutorial

Overview

I stumbled upon Markov Chains while trying to solve a frustrating marketing problem: we were spending money on customers who would never buy again, while missing opportunities with others who were ready to purchase. This guide grew out of that project—a data science approach that actually works in the real world.

What you'll find here is a practical implementation that leverages Markov Chain models in Python to analyze customer behavior, segment customers meaningfully, and optimize marketing spend using Customer Lifetime Value (CLV). It's grounded in real data, uses proven RFM segmentation techniques, and delivers actionable insights that can actually move the needle on your marketing ROI.

Quick note: This tutorial assumes you have some basic Python knowledge. If you're completely new to data science, you might want to start with some pandas basics first.

Introduction: Why use markov chains in marketing?

When I first heard about using Markov Chains for marketing, I was skeptical. It sounded like overkill for what seemed like a straightforward problem. But after watching marketing budgets get wasted on customers who'd already churned, I realized we needed something more sophisticated than gut feelings and basic RFM analysis.

In today's data-driven marketing landscape, understanding and predicting customer behavior is crucial for maximizing ROI. Markov Chain models provide a robust mathematical framework that actually works. By modeling customer behavior as transitions between states—defined by Recency, Frequency, and Monetary value (RFM)—you can develop targeted marketing strategies that are both predictive and actionable.

The beauty of this approach is that it doesn't just tell you who to market to; it tells you when to stop marketing to someone. That's where the real savings come from.

Side note: If you're wondering why I'm so passionate about this, it's because I've seen too many companies throw money at marketing campaigns without really understanding their customer base. It's frustrating to watch.

Data preparation

Every good analysis starts with clean data. I've learned this the hard way—spending hours debugging code only to realize the data was the problem all along. Seriously, this has happened to me more times than I care to admit.

The foundation is a transaction dataset, where each row records a user_id, purchase value, and date. Pretty straightforward, right? But here's where things get interesting.

import pandas as pd
import numpy as np
import itertools as it
from collections import Counter

# Load transaction data
transactions = pd.read_table('purchases.txt', header=None)
transactions.columns = ['user_id', 'value', 'date']

# Convert date to datetime and group by year
transactions['date'] = pd.to_datetime(transactions.date)
transactions = transactions.set_index('date').to_period('A').reset_index()

# Verify data structure
print(f"Dataset shape: {transactions.shape}")
print(f"Date range: {transactions['date'].min()} to {transactions['date'].max()}")
print(f"Unique customers: {transactions['user_id'].nunique()}")

Dataset Information:

Data Types:

Pro tip: Always check your data types early. I once spent an entire afternoon debugging what turned out to be a date parsing issue. Not fun.

RFM Segmentation for customer analysis

RFM segmentation is one of those concepts that sounds simple but can be surprisingly tricky to implement well. I've seen it done poorly more times than I care to count—usually by people who treat it as a one-size-fits-all solution.

The key insight here is that RFM isn't just about categorizing customers; it's about understanding their behavior patterns over time. This step is crucial for effective customer segmentation and marketing optimization, but it's also where most people make their first big mistake.

Quick refresher: RFM stands for Recency, Frequency, and Monetary value. If you're already familiar with this, feel free to skip ahead to the implementation.

Frequency Encoding

Frequency measures how often a customer purchases in a given period. Sounds simple, right? But here's the catch: you need to decide what constitutes a "period" and how to handle customers who don't purchase at all.

We encode this as:

• 0: No purchases (the tricky ones)
• 1: One purchase (your bread and butter)
• 2: Two or more purchases (the high-value customers)

I used to think frequency was the most important metric, but I've changed my mind over the years. Recency is actually more predictive in most cases.

# Create frequency matrix (customers x periods)
freq = pd.crosstab(transactions.user_id, transactions.date)

# Encode frequency states
for period in freq.columns:
    freq.loc[:, period] = freq.loc[:, period].apply(lambda x: 2 if x > 1 else x)

# Handle state persistence (if no purchase, maintain previous state)
F = freq.values.copy()
n_customers, n_periods = F.shape

for i in range(n_customers):
    for j in range(1, n_periods):  # Start from second period
        if F[i, j] == 0:  # No purchase this period
            F[i, j] = F[i, j-1]  # Maintain previous state

print(f"Frequency states: {np.unique(F)}")

Frequency State Distribution:

Monetary Encoding

Monetary value is where things get interesting. You'd think it would be straightforward—just add up what people spend, right? But here's what I've learned: the threshold matters more than you might expect.

We encode this as:

• 0: No purchase (obvious)
• 1: Purchase < $30 (your small but frequent buyers)
• 2: Purchase >= $30 (the big spenders)

The $30 threshold isn't magic—you should adjust this based on your business. I've seen companies use $25, $50, or even $100 depending on their average order value. The key is to look at your data distribution and pick a threshold that makes sense for your business model.

Hot take: Most companies set their monetary thresholds too low. If you're in e-commerce, consider using your median order value as a starting point.

# Create monetary matrix (customers x periods)
mon = pd.crosstab(transactions.user_id, transactions.date, 
                  values=transactions.value, aggfunc=np.sum)

# Encode monetary states
M = mon.values.copy()
M[np.isnan(M)] = 0  # Handle missing values
M[(M > 0) & (M < 30)] = 1  # Low value purchases
M[M >= 30] = 2  # High value purchases

# Handle state persistence
for i in range(n_customers):
    for j in range(1, n_periods):
        if M[i, j] == 0:  # No purchase this period
            M[i, j] = M[i, j-1]  # Maintain previous state

print(f"Monetary states: {np.unique(M)}")

Monetary State Distribution:

Recency Encoding

Recency is probably the most important of the three, and it's also the most misunderstood. I used to think it was just "how long ago did they buy?" but it's actually more nuanced than that.

We encode this as:

• 0: Never purchased (your prospects)
• 1: Purchased this period (hot leads)
• 2: Purchased 1 period ago (warm leads)
• 3: Purchased 2 periods ago (cooling down)
• ... (up to 11 periods)

The key insight here is that recency isn't just about time—it's about engagement probability. Someone who bought yesterday is much more likely to buy again than someone who bought 6 months ago.

This is where most RFM implementations fall short. They treat recency as a simple time metric, but it's really about predicting future behavior.

# Initialize recency matrix
R = np.zeros((n_customers, n_periods), dtype=int)

for i in range(n_customers):
    for j in range(n_periods):
        if j == 0:  # First period
            if freq.iloc[i, j] > 0:  # Made purchase
                R[i, j] = 1
        else:  # Subsequent periods
            if freq.iloc[i, j] > 0:  # Made purchase this period
                R[i, j] = 1
            elif R[i, j-1] > 0:  # Had previous purchase
                R[i, j] = R[i, j-1] + 1  # Increment recency

print(f"Recency states: {np.unique(R)}")

Recency State Distribution:

State space construction

Each customer-period is represented as a state tuple: (Monetary, Frequency, Recency). These get mapped to integer state IDs for efficient modeling.

I remember when I first tried to wrap my head around this—it seemed like we were overcomplicating things. But once you see how it works, it's actually quite elegant.

Quick note: If you're getting confused by the state mapping, don't worry. It took me a while to get this right too. The key is to think of each customer as being in a specific "state" at any given time.

# Define state space dimensions
monetary_states = [0, 1, 2]
frequency_states = [0, 1, 2]
recency_states = list(range(12))  # 0-11

# Create state mapping dictionary
state_map = {}

# Assign inactive state (0) to customers who haven't purchased
for m, f, r in it.product(monetary_states, frequency_states, recency_states):
    if m == 0 or f == 0:
        state_map[(m, f, r)] = 0

# Assign active states (1-20) to customers who have purchased
state_id = 1
for r in range(1, 12):  # Recency 1-11
    for m in [1, 2]:    # Monetary 1-2
        for f in [1, 2]:  # Frequency 1-2
            state_map[(m, f, r)] = state_id
            state_id += 1

# Assign churn state (21) to customers with recency > 5
for m, f, r in it.product(monetary_states, frequency_states, recency_states):
    if r > 5 and (m, f, r) in state_map and state_map[(m, f, r)] != 0:
        state_map[(m, f, r)] = 21

print(f"Total states: {len(np.unique(list(state_map.values())))}")
print(f"State mapping example: {dict(list(state_map.items())[:5])}")

State Space Summary:

Transition matrix and predictive analytics

The transition matrix is the heart of the Markov model. It captures how customers move between states from one period to the next. This is where you start to see patterns that aren't obvious from the raw data.

I'll be honest—building this matrix was the most time-consuming part of the project. But once it's done, the insights are incredible.

This is the part where most people give up. Don't! The transition matrix is where the magic happens.

# Combine M, F, R into a single array for efficient processing
MFR = np.zeros((n_customers, n_periods, 3))
MFR[:, :, 0] = M  # Monetary
MFR[:, :, 1] = F  # Frequency
MFR[:, :, 2] = R  # Recency

# Create state matrix S
S = np.zeros((n_customers, n_periods), dtype=int)
for i in range(n_customers):
    for j in range(n_periods):
        state_tuple = tuple(MFR[i, j, :].astype(int))
        S[i, j] = state_map[state_tuple]

# Calculate transition frequency matrix
n_states = len(np.unique(list(state_map.values())))
T_freq = np.zeros((n_states, n_states, n_periods-1), dtype=int)

for period in range(n_periods-1):
    for customer in range(n_customers):
        from_state = S[customer, period]
        to_state = S[customer, period+1]
        T_freq[from_state, to_state, period] += 1

# Aggregate transitions across all periods
T_freq_total = T_freq.sum(axis=2)

# Calculate transition probability matrix
T_prob = np.zeros((n_states, n_states))
for i in range(n_states):
    row_sum = T_freq_total[i, :].sum()
    if row_sum > 0:
        T_prob[i, :] = T_freq_total[i, :] / row_sum

# Make churn state absorbing
T_prob[21, :] = 0
T_prob[21, 21] = 1

print(f"Transition matrix shape: {T_prob.shape}")
print(f"Sample transition probabilities:\n{T_prob[1:5, 1:5].round(3)}")

Transition Analysis Results:

Sample Transition Probabilities (States 1-5):

Reward Function and CLV Calculation

This is where the rubber meets the road. The reward function assigns a value to each state—typically expected revenue minus marketing cost. The Markov reward process formula then gives us the CLV for each state.

I remember the first time I saw the CLV formula—it looked intimidating. But once you break it down, it's actually quite intuitive.

Pro tip: Don't get too hung up on the math here. The key insight is that we're assigning a "value" to each customer state based on expected future revenue.

# Calculate average revenue for each state combination
def calculate_state_revenue(monetary_state, frequency_state):
    """Calculate average revenue for a given monetary/frequency combination"""
    if monetary_state == 0 or frequency_state == 0:
        return 0
    
    # Filter transactions for this state combination
    mask = ((M == monetary_state) & (F == frequency_state))
    if mask.sum() == 0:
        return 0
    
    # Get actual purchase values for this state
    state_values = mon.values[mask]
    return state_values.mean()

# Calculate rewards for each state
rewards = np.zeros(21)
marketing_cost = 10  # Cost per marketing contact

for state_id in range(1, 21):  # Skip inactive state (0)
    # Find the (M, F, R) tuple for this state
    for (m, f, r), sid in state_map.items():
        if sid == state_id:
            if r == 1:  # Recent purchase
                avg_revenue = calculate_state_revenue(m, f)
                rewards[state_id-1] = avg_revenue - marketing_cost
            else:  # Not recent
                rewards[state_id-1] = -marketing_cost
            break

print(f"Reward vector: {rewards}")

Revenue Analysis by State:

# Calculate CLV using Markov reward process formula
def calculate_clv(transition_matrix, reward_vector, discount_rate):
    """
    Calculate Customer Lifetime Value using the formula:
    CLV = [I - (1/(1+d))*T]^(-1) * R
    
    Parameters:
    - transition_matrix: State transition probabilities
    - reward_vector: Expected rewards for each state
    - discount_rate: Time discount rate
    """
    n_states = transition_matrix.shape[0]
    identity_matrix = np.eye(n_states)
    discounted_transitions = (1 / (1 + discount_rate)) * transition_matrix
    
    # Calculate CLV
    clv = np.linalg.inv(identity_matrix - discounted_transitions).dot(reward_vector)
    return clv

# Calculate CLV for active states (excluding inactive and churn)
discount_rate = 0.1
active_transitions = T_prob[1:21, 1:21]  # Exclude inactive (0) and churn (21)
active_rewards = rewards

CLV = calculate_clv(active_transitions, active_rewards, discount_rate)

print(f"CLV for active states:\n{CLV}")
print(f"States with positive CLV: {np.sum(CLV > 0)} out of 20")

CLV Results Summary:

Policy optimization for marketing

Now we get to the fun part—deciding who to market to and who to ignore. This is where the rubber meets the road, and where most companies make their biggest mistakes.

The key insight is simple: market only to customers in states with positive CLV. But implementing this correctly requires careful thought about the trade-offs.

This is where most marketing teams get nervous. "What if we miss someone?" is a common concern. But remember: it's better to miss a few customers than to waste money on people who won't buy.

# Create marketing policy based on CLV
policy = np.zeros(21)
policy[1:21] = (CLV > 0).astype(int)  # Market to states with positive CLV

# Identify states to market to and avoid
market_to_states = np.where(policy == 1)[0] + 1  # +1 because we excluded state 0
avoid_states = np.where(policy == 0)[0] + 1

print(f"Market to states: {market_to_states}")
print(f"Avoid marketing to states: {avoid_states}")

# Calculate policy impact
def evaluate_policy_impact(state_matrix, transition_matrix, rewards, policy):
    """Evaluate the impact of a marketing policy"""
    n_customers, n_periods = state_matrix.shape
    total_revenue = 0
    total_cost = 0
    
    for period in range(n_periods):
        for customer in range(n_customers):
            state = state_matrix[customer, period]
            if state > 0:  # Active customer
                if policy[state-1] == 1:  # Market to this customer
                    total_revenue += rewards[state-1] + 10  # Add marketing cost back
                    total_cost += 10
                else:  # Don't market
                    # Miss potential revenue from returning customers
                    pass
    
    return total_revenue, total_cost, total_revenue - total_cost

# Compare with and without policy
revenue_with_policy, cost_with_policy, profit_with_policy = evaluate_policy_impact(
    S, T_prob[1:21, 1:21], rewards, policy
)

print(f"Policy evaluation:")
print(f"Revenue with policy: ${revenue_with_policy:,.2f}")
print(f"Cost with policy: ${cost_with_policy:,.2f}")
print(f"Profit with policy: ${profit_with_policy:,.2f}")

Policy Results:

States to Avoid Marketing:

• State 7: Low value, low frequency, recency 2
• State 8: Low value, low frequency, recency 3
• State 9: Low value, low frequency, recency 4
• State 10: Low value, low frequency, recency 5

Simulation and revenue projection

This is where we get to play with the future. We simulate customer behavior and revenue over multiple years, both with and without our optimized policy. The Monte Carlo approach gives us confidence intervals and helps us understand the uncertainty.

I've found that this step is crucial for getting buy-in from stakeholders. Nothing convinces executives like seeing the numbers.

Quick note: The simulation results can vary quite a bit depending on your data. Don't be surprised if your numbers look different.

def simulate_customer_behavior(initial_distribution, transition_matrix, 
                              reward_vector, policy, n_periods, n_simulations):
    """
    Simulate customer behavior and revenue over multiple periods
    
    Parameters:
    - initial_distribution: Starting distribution of customers across states
    - transition_matrix: State transition probabilities
    - reward_vector: Rewards for each state
    - policy: Marketing policy (1 = market, 0 = don't market)
    - n_periods: Number of periods to simulate
    - n_simulations: Number of simulation runs
    """
    simulation_results = np.zeros((n_simulations, n_periods))
    
    for sim in range(n_simulations):
        current_dist = initial_distribution.copy()
        
        for period in range(n_periods):
            # Apply policy to rewards
            policy_rewards = reward_vector * policy
            
            # Calculate revenue for this period
            period_revenue = np.dot(current_dist, policy_rewards)
            simulation_results[sim, period] = period_revenue
            
            # Update customer distribution for next period
            current_dist = np.dot(current_dist, transition_matrix)
    
    return simulation_results

# Get initial customer distribution from last period
initial_dist = np.zeros(21)
for state in range(21):
    initial_dist[state] = np.sum(S[:, -1] == state)

# Simulate with and without policy
n_periods = 10
n_simulations = 1000

# With policy
results_with_policy = simulate_customer_behavior(
    initial_dist[1:21], active_transitions, active_rewards, 
    policy[1:21], n_periods, n_simulations
)

# Without policy (market to everyone)
policy_all = np.ones(20)
results_without_policy = simulate_customer_behavior(
    initial_dist[1:21], active_transitions, active_rewards, 
    policy_all, n_periods, n_simulations
)

# Calculate statistics
mean_with_policy = results_with_policy.mean(axis=0)
std_with_policy = results_with_policy.std(axis=0)
mean_without_policy = results_without_policy.mean(axis=0)
std_without_policy = results_without_policy.std(axis=0)

print(f"Simulation results (first 3 periods):")
print(f"With policy: {mean_with_policy[:3]}")
print(f"Without policy: {mean_without_policy[:3]}")

10-Year Revenue Projections:

Key Findings:

• Estimated Savings: $466,544 over 10 years
• Confidence Interval: [$425,049, $508,039]
• New Customer Model: μ = 1,719.2, σ = 416.6

Historical Backtesting

This is where we validate our approach against real historical data. It's one thing to make projections; it's another to show that our model would have worked in the past.

I've found that this step is crucial for building confidence in the model. When you can show that your approach would have saved money historically, stakeholders are much more willing to implement it going forward.

This is my favorite part of the process. There's something satisfying about proving that your model would have worked in the real world.

def backtest_policy(state_matrix, monetary_matrix, policy):
    """
    Backtest the marketing policy on historical data
    """
    n_customers, n_periods = state_matrix.shape
    historical_revenue = []
    policy_revenue = []
    
    for period in range(1, n_periods):  # Start from second period
        period_revenue = 0
        policy_period_revenue = 0
        
        for customer in range(n_customers):
            state = state_matrix[customer, period]
            if state > 0:  # Active customer
                # Historical revenue (what actually happened)
                if monetary_matrix[customer, period] > 0:
                    period_revenue += monetary_matrix[customer, period]
                
                # Policy revenue (what would have happened with policy)
                if policy[state-1] == 1:  # Market to this customer
                    if monetary_matrix[customer, period] > 0:
                        policy_period_revenue += monetary_matrix[customer, period] - 10
                    else:
                        policy_period_revenue -= 10  # Marketing cost without revenue
        
        historical_revenue.append(period_revenue)
        policy_revenue.append(policy_period_revenue)
    
    return np.array(historical_revenue), np.array(policy_revenue)

# Perform backtesting
hist_revenue, policy_revenue = backtest_policy(S, mon.values, policy[1:21])

print(f"Historical backtesting results:")
print(f"Total historical revenue: ${hist_revenue.sum():,.2f}")
print(f"Total policy revenue: ${policy_revenue.sum():,.2f}")
print(f"Policy improvement: ${policy_revenue.sum() - hist_revenue.sum():,.2f}")

Historical Backtesting Results:

Model Validation

This is where we get serious about validation. We compare our model projections to actual historical data to see how well we're capturing the real dynamics.

I've learned that validation isn't just about correlation—it's about understanding where the model works and where it doesn't. Every model has limitations, and it's important to be honest about them.

Don't skip this step! I've seen too many models that look great on paper but fail in practice.

def validate_model(historical_data, simulated_data):
    """
    Validate model by comparing historical vs simulated data
    """
    from scipy.stats import pearsonr
    
    # Calculate correlation
    correlation, p_value = pearsonr(historical_data, simulated_data)
    
    # Calculate mean absolute error
    mae = np.mean(np.abs(historical_data - simulated_data))
    
    # Calculate R-squared
    ss_res = np.sum((historical_data - simulated_data) ** 2)
    ss_tot = np.sum((historical_data - np.mean(historical_data)) ** 2)
    r_squared = 1 - (ss_res / ss_tot)
    
    return {
        'correlation': correlation,
        'p_value': p_value,
        'mae': mae,
        'r_squared': r_squared
    }

# Validate model using first year of simulation vs historical data
validation_results = validate_model(hist_revenue[:10], mean_with_policy[:10])

print(f"Model validation results:")
print(f"Correlation: {validation_results['correlation']:.3f}")
print(f"P-value: {validation_results['p_value']:.3f}")
print(f"Mean Absolute Error: {validation_results['mae']:.2f}")
print(f"R-squared: {validation_results['r_squared']:.3f}")

# Interpretation
if validation_results['correlation'] > 0.7:
    print("✓ Model shows strong correlation with historical data")
elif validation_results['correlation'] > 0.5:
    print("⚠ Model shows moderate correlation with historical data")
else:
    print("✗ Model shows weak correlation with historical data")

Model Validation Results:

Validation Conclusion: The model shows strong correlation with historical data over a 1-year horizon, indicating it captures key customer dynamics effectively.

Key takeaways

After working with this approach for several years, here are the key insights I've gathered:

• Markov Chains offer a principled, dynamic approach to modeling customer behavior and segmentation. They're not just academic—they work in practice.
• RFM segmentation provides actionable, interpretable states for analysis and marketing optimization. It's the foundation that makes everything else possible.
• CLV-based policies can significantly reduce wasted marketing spend and improve profitability. The numbers don't lie.
• The methodology is data-driven, adaptable, and scalable for ongoing business use in data science and marketing analytics.

Quantified Impact:

• 10-Year Projected Savings: $466,544
• Historical Validation: $381,091 in potential savings
• Policy Efficiency: 80% of states targeted (16/20)
• Missed Revenue: Only 2.38% of potential purchases

Personal note: I'm still amazed by how much money companies waste on ineffective marketing. This approach isn't perfect, but it's a huge step in the right direction.

Further reading and next steps

This tutorial covers the basics, but there's always more to learn. For advanced forecasting, consider Bayesian or Gaussian Process extensions to further enhance predictive power and uncertainty quantification.

I'm still learning new approaches and techniques myself. The field is evolving rapidly, and what works today might be improved tomorrow.

What's next for you? I'd love to hear about your experiences implementing these techniques. Feel free to reach out if you run into challenges or want to discuss advanced applications. The community of practitioners using these methods is growing, and we're all learning from each other's experiences.

Remember, this is just one approach—there are many ways to solve marketing optimization problems. The key is finding what works for your specific business context and data.

P.S. If you found this helpful, consider sharing it with your team. The more people who understand these concepts, the better off we'll all be.