1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

import numpy as np
from typing import List, Dict, Tuple
from dataclasses import dataclass
from scipy.optimize import minimize
 
@dataclass
class Impression:
    """Single impression with position and outcome."""
    item_id: str
    position: int
    clicked: bool
    user_id: str
 
 
class PositionBiasEstimator:
    """
    Estimate position bias from click data.
    
    Uses Expectation-Maximization to jointly estimate:
    - Position examination probabilities
    - Item relevance probabilities
    """
    
    def __init__(self, n_positions: int = 10):
        self.n_positions = n_positions
        self.examination_probs: np.ndarray = None
        self.relevance_scores: Dict[str, float] = {}
    
    def fit(
        self,
        impressions: List[Impression],
        max_iterations: int = 100,
        tol: float = 1e-4,
    ):
        """
        Estimate position bias using EM algorithm.
        
        Click probability = P(examine) × P(relevant)
        
        E-step: estimate click was due to examination given model
        M-step: update examination and relevance parameters
        """
        # Initialize
        self.examination_probs = np.array([
            1.0 / (i + 1) for i in range(self.n_positions)
        ])
        self.examination_probs /= self.examination_probs[0]  # Normalize
        
        # Initial relevance: empirical CTR per item
        item_clicks = {}
        item_impressions = {}
        for imp in impressions:
            item_clicks[imp.item_id] = item_clicks.get(imp.item_id, 0) + imp.clicked
            item_impressions[imp.item_id] = item_impressions.get(imp.item_id, 0) + 1
        
        for item_id in item_impressions:
            self.relevance_scores[item_id] = (
                item_clicks.get(item_id, 0) / item_impressions[item_id]
            )
        
        # EM iterations
        prev_ll = float('-inf')
        for iteration in range(max_iterations):
            # E-step: compute expected examination given click
            exam_expectations = self._e_step(impressions)
            
            # M-step: update parameters
            self._m_step(impressions, exam_expectations)
            
            # Check convergence
            ll = self._log_likelihood(impressions)
            if ll - prev_ll < tol:
                print(f"Converged at iteration {iteration}")
                break
            prev_ll = ll
    
    def _e_step(
        self,
        impressions: List[Impression],
    ) -> List[float]:
        """Compute expected examination probability."""
        expectations = []
        
        for imp in impressions:
            exam_prob = self.examination_probs[min(imp.position, self.n_positions - 1)]
            rel_prob = self.relevance_scores.get(imp.item_id, 0.5)
            
            if imp.clicked:
                # Click implies examination
                expectations.append(1.0)
            else:
                # No click: could be not examined or not relevant
                # P(not examined | not clicked) using Bayes
                p_not_click_not_exam = 1 - exam_prob
                p_not_click_exam_not_rel = exam_prob * (1 - rel_prob)
                p_not_click = p_not_click_not_exam + p_not_click_exam_not_rel
                
                p_exam_given_not_click = p_not_click_exam_not_rel / max(p_not_click, 1e-10)
                expectations.append(p_exam_given_not_click)
        
        return expectations
    
    def _m_step(
        self,
        impressions: List[Impression],
        exam_expectations: List[float],
    ):
        """Update examination and relevance parameters."""
        # Update examination probabilities
        position_exam_sum = np.zeros(self.n_positions)
        position_count = np.zeros(self.n_positions)
        
        for imp, exam_exp in zip(impressions, exam_expectations):
            pos = min(imp.position, self.n_positions - 1)
            position_exam_sum[pos] += exam_exp
            position_count[pos] += 1
        
        for pos in range(self.n_positions):
            if position_count[pos] > 0:
                self.examination_probs[pos] = position_exam_sum[pos] / position_count[pos]
        
        # Update relevance scores
        item_click_sum = {}
        item_exam_sum = {}
        
        for imp, exam_exp in zip(impressions, exam_expectations):
            item_click_sum[imp.item_id] = (
                item_click_sum.get(imp.item_id, 0) + imp.clicked
            )
            item_exam_sum[imp.item_id] = (
                item_exam_sum.get(imp.item_id, 0) + exam_exp
            )
        
        for item_id in item_exam_sum:
            if item_exam_sum[item_id] > 0:
                self.relevance_scores[item_id] = (
                    item_click_sum.get(item_id, 0) / item_exam_sum[item_id]
                )
    
    def _log_likelihood(self, impressions: List[Impression]) -> float:
        """Compute log likelihood of data given current parameters."""
        ll = 0.0
        for imp in impressions:
            pos = min(imp.position, self.n_positions - 1)
            exam_prob = self.examination_probs[pos]
            rel_prob = self.relevance_scores.get(imp.item_id, 0.5)
            
            click_prob = exam_prob * rel_prob
            
            if imp.clicked:
                ll += np.log(max(click_prob, 1e-10))
            else:
                ll += np.log(max(1 - click_prob, 1e-10))
        
        return ll
    
    def get_debiased_relevance(self, item_id: str) -> float:
        """Get position-debiased relevance estimate."""
        return self.relevance_scores.get(item_id, 0.5)
    
    def get_position_bias_curve(self) -> np.ndarray:
        """Get estimated examination probability by position."""
        return self.examination_probs.copy()
 
 
class IPSEstimator:
    """
    Inverse Propensity Scoring for selection bias correction.
    
    Reweights observations by inverse of selection probability
    to get unbiased estimates.
    """
    
    def __init__(self, propensity_model=None):
        """
        Args:
            propensity_model: Model that predicts P(item shown | context)
                             If None, uses historical frequencies.
        """
        self.propensity_model = propensity_model
        self.item_propensities: Dict[str, float] = {}
    
    def estimate_propensities_from_logs(
        self,
        impressions: List[Impression],
        total_opportunities: int = None,
    ):
        """
        Estimate propensities from historical impression logs.
        
        P(item shown) ≈ impressions / total_opportunities
        """
        item_counts = {}
        for imp in impressions:
            item_counts[imp.item_id] = item_counts.get(imp.item_id, 0) + 1
        
        if total_opportunities is None:
            total_opportunities = len(impressions)
        
        for item_id, count in item_counts.items():
            self.item_propensities[item_id] = count / total_opportunities
    
    def compute_ips_estimate(
        self,
        impressions: List[Impression],
        clip_threshold: float = 0.01,
    ) -> Dict[str, float]:
        """
        Compute IPS-weighted relevance estimates.
        
        IPS weights: w_i = 1 / P(shown_i)
        
        With clipping to avoid extreme weights for rare items.
        """
        item_weighted_clicks = {}
        item_weights = {}
        
        for imp in impressions:
            propensity = max(
                self.item_propensities.get(imp.item_id, 0.5),
                clip_threshold  # Clip to avoid huge weights
            )
            weight = 1.0 / propensity
            
            item_weighted_clicks[imp.item_id] = (
                item_weighted_clicks.get(imp.item_id, 0) + 
                weight * imp.clicked
            )
            item_weights[imp.item_id] = (
                item_weights.get(imp.item_id, 0) + weight
            )
        
        # Normalize
        estimates = {}
        for item_id in item_weights:
            estimates[item_id] = (
                item_weighted_clicks[item_id] / item_weights[item_id]
            )
        
        return estimates

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

import numpy as np
from typing import List, Dict, Tuple, Callable
from dataclasses import dataclass
 
@dataclass
class LoggedInteraction:
    """Single logged interaction from production."""
    context: np.ndarray  # User/session features
    action: str  # Item shown
    reward: float  # Click/conversion
    propensity: float  # P(action | context) under logging policy
 
 
class CounterfactualEvaluator:
    """
    Evaluate new recommendation policies using logged data.
    
    Implements IPS, SNIPS, and Doubly Robust estimators.
    """
    
    def __init__(
        self,
        weight_clip: float = 10.0,
        min_propensity: float = 0.01,
    ):
        self.weight_clip = weight_clip
        self.min_propensity = min_propensity
    
    def evaluate_ips(
        self,
        logged_data: List[LoggedInteraction],
        target_policy: Callable[[np.ndarray, str], float],
    ) -> Tuple[float, float]:
        """
        Inverse Propensity Scoring estimation.
        
        Args:
            logged_data: Interactions from logging policy
            target_policy: Function(context, action) -> P(action | context)
                          under the target policy
        
        Returns:
            (estimate, standard_error)
        """
        weighted_rewards = []
        
        for interaction in logged_data:
            # Propensity ratio
            target_prob = target_policy(interaction.context, interaction.action)
            logging_prob = max(interaction.propensity, self.min_propensity)
            
            weight = target_prob / logging_prob
            weight = min(weight, self.weight_clip)  # Clip
            
            weighted_rewards.append(weight * interaction.reward)
        
        estimate = np.mean(weighted_rewards)
        std_error = np.std(weighted_rewards) / np.sqrt(len(weighted_rewards))
        
        return estimate, std_error
    
    def evaluate_snips(
        self,
        logged_data: List[LoggedInteraction],
        target_policy: Callable[[np.ndarray, str], float],
    ) -> Tuple[float, float]:
        """
        Self-Normalized IPS estimation.
        
        Normalizes by sum of weights for lower variance.
        """
        weights = []
        weighted_rewards = []
        
        for interaction in logged_data:
            target_prob = target_policy(interaction.context, interaction.action)
            logging_prob = max(interaction.propensity, self.min_propensity)
            
            weight = min(target_prob / logging_prob, self.weight_clip)
            weights.append(weight)
            weighted_rewards.append(weight * interaction.reward)
        
        total_weight = sum(weights)
        estimate = sum(weighted_rewards) / total_weight if total_weight > 0 else 0
        
        # Approximate standard error
        n = len(logged_data)
        variance = np.var([
            w * (r - estimate) for w, r 
            in zip(weights, [i.reward for i in logged_data])
        ]) / (total_weight ** 2)
        std_error = np.sqrt(variance * n)
        
        return estimate, std_error
    
    def evaluate_doubly_robust(
        self,
        logged_data: List[LoggedInteraction],
        target_policy: Callable[[np.ndarray, str], float],
        reward_model: Callable[[np.ndarray, str], float],
    ) -> Tuple[float, float]:
        """
        Doubly Robust estimation.
        
        Combines direct method with IPS correction.
        Consistent if either reward_model or propensities are correct.
        
        Args:
            reward_model: Function(context, action) -> E[reward]
        """
        estimates = []
        
        for interaction in logged_data:
            target_prob = target_policy(interaction.context, interaction.action)
            logging_prob = max(interaction.propensity, self.min_propensity)
            
            weight = min(target_prob / logging_prob, self.weight_clip)
            
            # Predicted reward for logged action
            predicted_reward = reward_model(interaction.context, interaction.action)
            
            # Expected reward under target policy (approximate)
            # In practice, would sum over all actions weighted by target_prob
            expected_under_target = reward_model(
                interaction.context, 
                interaction.action  # Simplified: same action
            )
            
            # DR estimate
            dr_term = expected_under_target + weight * (interaction.reward - predicted_reward)
            estimates.append(dr_term)
        
        estimate = np.mean(estimates)
        std_error = np.std(estimates) / np.sqrt(len(estimates))
        
        return estimate, std_error
 
 
class PolicyComparisonFramework:
    """
    Framework for comparing multiple recommendation policies
    using counterfactual evaluation.
    """
    
    def __init__(self):
        self.evaluator = CounterfactualEvaluator()
    
    def compare_policies(
        self,
        logged_data: List[LoggedInteraction],
        policies: Dict[str, Callable],
        methods: List[str] = ['ips', 'snips'],
    ) -> Dict[str, Dict[str, Tuple[float, float]]]:
        """
        Compare multiple policies using multiple estimators.
        
        Returns:
            {policy_name: {method: (estimate, stderr)}}
        """
        results = {}
        
        for policy_name, policy_fn in policies.items():
            results[policy_name] = {}
            
            for method in methods:
                if method == 'ips':
                    est, se = self.evaluator.evaluate_ips(logged_data, policy_fn)
                elif method == 'snips':
                    est, se = self.evaluator.evaluate_snips(logged_data, policy_fn)
                else:
                    continue
                
                results[policy_name][method] = (est, se)
                print(f"{policy_name} ({method}): {est:.4f} ± {se:.4f}")
        
        return results
    
    def compute_confidence_interval(
        self,
        estimate: float,
        std_error: float,
        confidence: float = 0.95,
    ) -> Tuple[float, float]:
        """Compute confidence interval for estimate."""
        from scipy import stats
        z = stats.norm.ppf((1 + confidence) / 2)
        return (estimate - z * std_error, estimate + z * std_error)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

import numpy as np
from typing import List, Dict, Tuple
from collections import defaultdict
 
class FeedbackLoopMonitor:
    """
    Monitor for detecting feedback loop symptoms.
    
    Tracks metrics over time to identify concerning trends.
    """
    
    def __init__(self, window_size: int = 7):
        self.window_size = window_size
        self.metric_history: Dict[str, List[float]] = defaultdict(list)
    
    def log_metrics(self, metrics: Dict[str, float], timestamp: str):
        """Log daily/periodic metrics for trend analysis."""
        for name, value in metrics.items():
            self.metric_history[name].append(value)
            # Keep only recent window
            if len(self.metric_history[name]) > self.window_size * 4:
                self.metric_history[name] = self.metric_history[name][-self.window_size * 4:]
    
    def detect_concentration(self) -> Dict[str, float]:
        """
        Detect increasing concentration (decreasing diversity).
        
        Returns trend coefficients for concentration metrics.
        """
        alerts = {}
        
        for metric in ['gini_coefficient', 'top_10_concentration', 'catalog_coverage']:
            if metric in self.metric_history and len(self.metric_history[metric]) >= self.window_size:
                values = self.metric_history[metric][-self.window_size:]
                
                # Compute trend
                x = np.arange(len(values))
                slope, _ = np.polyfit(x, values, 1)
                
                # Alert thresholds
                if metric == 'catalog_coverage' and slope < -0.01:  # Decreasing coverage
                    alerts[metric] = {'trend': slope, 'alert': 'Coverage decreasing'}
                elif metric == 'gini_coefficient' and slope > 0.01:  # Increasing inequality
                    alerts[metric] = {'trend': slope, 'alert': 'Concentration increasing'}
                elif metric == 'top_10_concentration' and slope > 0.005:
                    alerts[metric] = {'trend': slope, 'alert': 'Winner-take-all emerging'}
        
        return alerts
    
    def detect_filter_bubble(
        self,
        user_category_history: Dict[str, List[str]],
    ) -> Dict[str, float]:
        """
        Detect filter bubble formation.
        
        Measures user-level category diversity over time.
        """
        bubble_scores = {}
        
        for user_id, categories in user_category_history.items():
            if len(categories) < 10:
                continue
            
            # Split into early and recent
            early = categories[:len(categories)//2]
            recent = categories[len(categories)//2:]
            
            # Measure diversity (entropy)
            def entropy(items):
                if not items:
                    return 0
                unique, counts = np.unique(items, return_counts=True)
                probs = counts / len(items)
                return -np.sum(probs * np.log2(probs + 1e-10))
            
            early_entropy = entropy(early)
            recent_entropy = entropy(recent)
            
            # Declining entropy suggests bubble formation
            if early_entropy > 0:
                bubble_scores[user_id] = (early_entropy - recent_entropy) / early_entropy
        
        return bubble_scores
 
 
class LoopMitigationReranker:
    """
    Reranker that applies loop mitigation strategies.
    """
    
    def __init__(
        self,
        diversity_boost: float = 0.1,
        freshness_boost: float = 0.2,
        min_provider_slots: int = 3,
    ):
        self.diversity_boost = diversity_boost
        self.freshness_boost = freshness_boost
        self.min_provider_slots = min_provider_slots
        
        # Track impressions to boost underexposed items
        self.item_impressions: Dict[str, int] = defaultdict(int)
        self.total_impressions = 0
    
    def rerank(
        self,
        items: List[Dict],
        k: int,
    ) -> List[Dict]:
        """
        Rerank items with loop mitigation boosts.
        
        Items should have: item_id, score, provider_id, age_days
        """
        if not items:
            return []
        
        # Compute mitigation boosts
        for item in items:
            item['original_score'] = item['score']
            
            # Underexposure boost
            expected_impressions = self.total_impressions / max(len(self.item_impressions), 1)
            actual_impressions = self.item_impressions.get(item['item_id'], 0)
            
            if actual_impressions < expected_impressions * 0.5:
                underexposure_boost = self.diversity_boost
            else:
                underexposure_boost = 0
            
            # Freshness boost for new items
            age_days = item.get('age_days', 30)
            if age_days < 7:
                freshness_boost = self.freshness_boost * (1 - age_days / 7)
            else:
                freshness_boost = 0
            
            item['score'] = (
                item['original_score'] + 
                underexposure_boost + 
                freshness_boost
            )
        
        # Sort by boosted score
        items.sort(key=lambda x: x['score'], reverse=True)
        
        # Ensure provider diversity
        selected = []
        provider_counts = defaultdict(int)
        remaining = list(items)
        
        while len(selected) < k and remaining:
            for item in remaining:
                provider = item.get('provider_id')
                
                # Enforce max per provider
                if provider_counts[provider] < self.min_provider_slots:
                    selected.append(item)
                    remaining.remove(item)
                    provider_counts[provider] += 1
                    break
            else:
                # All providers at limit, just take next best
                if remaining:
                    selected.append(remaining.pop(0))
        
        return selected
    
    def record_impressions(self, item_ids: List[str]):
        """Record that items were shown."""
        for item_id in item_ids:
            self.item_impressions[item_id] += 1
            self.total_impressions += 1
 
 
class HistoricalDiversityMaintainer:
    """
    Maintain diversity by ensuring user sees items from
    categories they haven't seen recently.
    """
    
    def __init__(self, history_window: int = 50):
        self.history_window = history_window
        self.user_history: Dict[str, List[str]] = defaultdict(list)
    
    def boost_novel_categories(
        self,
        user_id: str,
        items: List[Dict],
    ) -> List[Dict]:
        """Boost items from categories user hasn't seen recently."""
        recent_categories = set(self.user_history[user_id][-self.history_window:])
        
        for item in items:
            category = item.get('category')
            if category and category not in recent_categories:
                item['score'] = item.get('score', 0) * 1.2  # 20% boost
        
        return items
    
    def record_shown(self, user_id: str, category: str):
        """Record shown category for user."""
        history = self.user_history[user_id]
        history.append(category)
        
        # Trim to window
        if len(history) > self.history_window * 2:
            self.user_history[user_id] = history[-self.history_window:]