Designing Sports Betting Systems in R: Bayesian Probabilities, Expected Value, and Kelly Logic

Posted on February 1, 2026 by rprogrammingbooks in R bloggers | 0 Comments

[This article was first published on Blog - R Programming Books, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A good sports betting system is not a “pick-winners” machine. It’s an uncertainty engine: it turns data into probabilities, probabilities into expected value, and expected value into position sizes that survive variance. If you can do those three steps consistently, you can build a robust process— even if individual bets lose often.

This post is a code-heavy blueprint for designing sports betting systems in R using Bayesian probabilities, market odds, EV, the Kelly Criterion, and practical betting strategies. You’ll also get a backtesting and Monte Carlo framework to evaluate drawdowns, ruin risk, and performance stability.

1. Setup: Packages, Data Schemas, and Conventions

We’ll work with an event-level dataset (one row per match) and optionally a bet-level dataset (one row per wager). The minimum fields you want for match win modeling:

date
team, opponent
home (1 home, 0 away)
result (1 win, 0 loss)
odds_decimal (decimal odds for the bet you want to evaluate)

# Core data + modeling stack
library(dplyr)
library(tidyr)
library(purrr)
library(stringr)
library(lubridate)
library(ggplot2)
library(scales)

# Bayesian regression (optional, heavier but powerful)
library(rstanarm)
library(posterior)

set.seed(42)

# --- Odds helpers ---
american_to_decimal <- function(american) {
  ifelse(american > 0, 1 + american/100, 1 + 100/abs(american))
}

decimal_to_american <- function(decimal) {
  # returns American odds as numeric (positive/negative)
  stopifnot(all(decimal > 1))
  out <- ifelse(decimal >= 2, (decimal - 1) * 100, -100/(decimal - 1))
  out
}

decimal_to_implied <- function(decimal) 1 / decimal

# Remove vig from a two-outcome market (simple proportional normalization)
remove_vig_two_way <- function(p1_raw, p2_raw) {
  s <- p1_raw + p2_raw
  c(p1 = p1_raw/s, p2 = p2_raw/s, overround = s)
}

# --- EV helpers ---
ev_decimal <- function(p, d) {
  # expected profit per unit stake for decimal odds d
  # win: profit = d - 1; lose: -1
  p*(d - 1) - (1 - p)
}

breakeven_p_decimal <- function(d) 1 / d

# --- Kelly helpers ---
kelly_fraction_decimal <- function(p, d) {
  # Kelly fraction for decimal odds d, win prob p
  # b = d - 1, q = 1-p, f* = (bp - q)/b
  b <- d - 1
  q <- 1 - p
  f <- (b*p - q) / b
  f
}

# Practical caps
clip <- function(x, lo, hi) pmax(lo, pmin(hi, x))

1.1 A toy dataset generator (so you can run everything)

In real life, you’ll load your match history and odds. For a self-contained tutorial, we simulate: latent team strengths → win probabilities → outcomes → bookmaker odds with a built-in margin.

simulate_league <- function(
  n_teams = 12,
  n_games = 2000,
  home_adv = 0.25,
  sigma_strength = 0.9,
  bookmaker_margin = 0.05
) {
  teams <- paste0("T", sprintf("%02d", 1:n_teams))
  strength <- rnorm(n_teams, 0, sigma_strength)
  names(strength) <- teams

  # schedule
  team <- sample(teams, n_games, replace = TRUE)
  opponent <- sample(teams, n_games, replace = TRUE)
  while(any(team == opponent)) {
    idx <- which(team == opponent)
    opponent[idx] <- sample(teams, length(idx), replace = TRUE)
  }

  home <- rbinom(n_games, 1, 0.5)
  date <- as.Date("2022-01-01") + sample(0:900, n_games, replace = TRUE)

  # true probability via logistic link
  lin <- (strength[team] - strength[opponent]) + home_adv * home
  p_true <- plogis(lin)

  # outcomes
  result <- rbinom(n_games, 1, p_true)

  # bookmaker sets odds from "market probability" with margin
  # market_prob ~ p_true with noise, then add margin by inflating implied probs
  p_mkt <- clip(p_true + rnorm(n_games, 0, 0.03), 0.02, 0.98)
  # apply margin by scaling implied probabilities upward (approx.)
  p_raw <- clip(p_mkt * (1 + bookmaker_margin), 0.02, 0.995)

  odds_decimal <- 1 / p_raw

  tibble(
    date = date,
    team = team,
    opponent = opponent,
    home = home,
    result = result,
    p_true = p_true,
    odds_decimal = odds_decimal
  ) %>%
    arrange(date)
}

matches <- simulate_league()
glimpse(matches)
summary(matches$odds_decimal)

2. Odds and Probability: Conversions, Vig (Overround), and Fair Prices

Odds imply probabilities, but sportsbook odds typically include a margin (vig/overround). If you compare your model’s probability to raw implied probability, you might misjudge edges. The right comparison is: your probability vs the market’s vig-free probability.

2.1 Basic conversions

# Example: decimal odds
d <- 1.91
p_implied_raw <- decimal_to_implied(d)
p_implied_raw
# Example: American odds
am <- -110
d2 <- american_to_decimal(am)
c(decimal = d2, implied = decimal_to_implied(d2))

2.2 Two-way vig removal example

For a two-outcome market (like moneyline without draws), you can normalize implied probabilities. Suppose the book offers:

d_home <- 1.80
d_away <- 2.10

p_home_raw <- 1/d_home
p_away_raw <- 1/d_away

vigfree <- remove_vig_two_way(p_home_raw, p_away_raw)
vigfree
# Market "fair" probabilities
p_home_fair <- vigfree["p1"]
p_away_fair <- vigfree["p2"]

p_home_fair + p_away_fair  # should be 1

In practice, you might not always have both sides’ odds, especially if you store only “the odds you bet”. But for sharp evaluation and backtesting, capturing the full market quote (both sides) is ideal.

3. Expected Value: Break-even Probability, Edge, and Sensitivity

Expected value (EV) is the bridge between probability and betting decisions. With decimal odds d, and win probability p, EV per 1 unit stake is:

# EV for a single bet
p_hat <- 0.56
d <- 1.95

ev <- ev_decimal(p_hat, d)
ev
# Break-even p (the book's implied probability)
p_be <- breakeven_p_decimal(d)
c(breakeven_p = p_be, edge = p_hat - p_be)

3.1 EV sensitivity curve

ev_curve <- function(d, p_grid = seq(0.01, 0.99, by = 0.01)) {
  tibble(p = p_grid, ev = ev_decimal(p_grid, d))
}

df_ev <- ev_curve(1.95)

ggplot(df_ev, aes(p, ev)) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  labs(
    title = "EV per Unit Stake vs Win Probability",
    x = "Win probability p",
    y = "Expected value (profit per 1 unit stake)"
  )

Notice how small probability errors can flip EV from positive to negative when odds are tight. This is why calibration and uncertainty-aware betting matter.

4. Bayesian Foundations for Betting: Priors, Posteriors, and Uncertainty

Bayesian thinking fits betting naturally: you start with a prior belief about a team/player/strategy, update with evidence, and end up with a posterior distribution over the probability you care about. Crucially, you get uncertainty, not just a point estimate.

4.1 A simple Bayesian coin: Beta prior + Binomial likelihood

# Prior: Beta(a, b)
a0 <- 10
b0 <- 10

# Observed: k wins out of n
k <- 12
n <- 20

# Posterior: Beta(a0 + k, b0 + n - k)
a1 <- a0 + k
b1 <- b0 + (n - k)

# Posterior mean and interval
post_mean <- a1 / (a1 + b1)
post_ci <- qbeta(c(0.05, 0.95), a1, b1)

c(mean = post_mean, lo90 = post_ci[1], hi90 = post_ci[2])
# Visualize prior vs posterior
grid <- seq(0, 1, length.out = 400)

df_beta <- bind_rows(
  tibble(p = grid, density = dbeta(grid, a0, b0), dist = "prior"),
  tibble(p = grid, density = dbeta(grid, a1, b1), dist = "posterior")
)

ggplot(df_beta, aes(p, density, color = dist)) +
  geom_line() +
  labs(
    title = "Beta Prior vs Posterior",
    x = "Win probability p",
    y = "Density"
  )

This is the simplest but most important idea: if your sample is small, your posterior won’t be overly confident. That alone can prevent overbetting.

5. Beta-Binomial Models: Bayesian Win Rates (Conjugate Bayes)

The Beta-Binomial is a workhorse for modeling win rates, conversion rates, “cover rates”, etc. You can use it for:

Team win rate vs certain opponents
Over/Under hit rate for a model-defined filter
Any binary market outcome you can label as 0/1

5.1 Team-level posteriors

# Build per-team counts
team_stats <- matches %>%
  group_by(team) %>%
  summarise(
    n = n(),
    k = sum(result),
    .groups = "drop"
  )

# Choose a prior: mildly skeptical
a0 <- 8
b0 <- 8

team_post <- team_stats %>%
  mutate(
    a = a0 + k,
    b = b0 + (n - k),
    post_mean = a/(a+b),
    lo90 = qbeta(0.05, a, b),
    hi90 = qbeta(0.95, a, b)
  ) %>%
  arrange(desc(post_mean))

team_post %>% slice(1:10)
ggplot(team_post, aes(reorder(team, post_mean), post_mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = lo90, ymax = hi90), width = 0.2) +
  coord_flip() +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  labs(
    title = "Team Win Probability Posteriors (Beta-Binomial)",
    x = "Team",
    y = "Posterior mean (90% interval)"
  )

5.2 Betting with uncertainty (posterior sampling)

Instead of using a single probability, sample from the posterior and compute EV distribution. This helps you avoid betting when your edge is too uncertain.

# Suppose you want to bet Team X at odds d
team_x <- team_post$team[1]
d <- 1.95

a <- team_post$a[team_post$team == team_x]
b <- team_post$b[team_post$team == team_x]

p_draws <- rbeta(20000, a, b)
ev_draws <- ev_decimal(p_draws, d)

quantile(ev_draws, c(0.05, 0.5, 0.95))
# Probability that EV is positive (posterior probability of an edge)
mean(ev_draws > 0)

A simple rule: only bet if P(EV > 0) exceeds a threshold (e.g., 0.6 or 0.7). That’s a Bayesian risk control.

6. Sequential Updating: Learning Over Time (Online Bayesian Updating)

Betting systems are online systems. Your beliefs should update as games happen. Beta-Binomial makes this trivial.

update_beta <- function(a, b, y) {
  # y is 1 for win, 0 for loss
  c(a = a + y, b = b + (1 - y))
}

sequential_posterior <- function(y, a0 = 8, b0 = 8) {
  a <- a0; b <- b0
  out <- vector("list", length(y))
  for (i in seq_along(y)) {
    ab <- update_beta(a, b, y[i])
    a <- ab["a"]; b <- ab["b"]
    out[[i]] <- c(
      i = i,
      a = a, b = b,
      mean = a/(a+b),
      lo90 = qbeta(0.05, a, b),
      hi90 = qbeta(0.95, a, b)
    )
  }
  bind_rows(lapply(out, as_tibble_row))
}

# Pick one team and track posterior over time
team_pick <- sample(unique(matches$team), 1)
y_seq <- matches %>% filter(team == team_pick) %>% arrange(date) %>% pull(result)

seq_df <- sequential_posterior(y_seq, a0 = 8, b0 = 8) %>%
  mutate(team = team_pick)

ggplot(seq_df, aes(i, mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lo90, ymax = hi90), alpha = 0.2) +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  labs(
    title = paste("Sequential Posterior for", team_pick),
    x = "Game index",
    y = "Posterior mean (90% band)"
  )

7. Shrinkage and Partial Pooling: Stabilizing Team/Player Estimates

Raw win rates are noisy. Fully separate Beta posteriors help, but you can go further: estimate the prior from the league (Empirical Bayes), then use it as a shrinkage prior for each team.

7.1 Empirical Bayes prior fit (method-of-moments)

fit_beta_mom <- function(p_hat, w = NULL) {
  # method of moments for Beta parameters
  # p_hat: observed proportions, w optional weights (e.g., n games)
  if (is.null(w)) w <- rep(1, length(p_hat))

  m <- weighted.mean(p_hat, w)
  v <- weighted.mean((p_hat - m)^2, w)

  # clamp variance to avoid weirdness
  v <- max(v, 1e-6)

  # Beta variance: m(1-m)/(a+b+1)
  t <- m*(1-m)/v - 1
  a <- max(0.5, m*t)
  b <- max(0.5, (1-m)*t)
  c(a = a, b = b)
}

p_hat <- team_stats$k / team_stats$n
prior_ab <- fit_beta_mom(p_hat, w = team_stats$n)
prior_ab
a0 <- prior_ab["a"]
b0 <- prior_ab["b"]

team_post_eb <- team_stats %>%
  mutate(
    a = a0 + k,
    b = b0 + (n - k),
    post_mean = a/(a+b),
    lo90 = qbeta(0.05, a, b),
    hi90 = qbeta(0.95, a, b)
  ) %>%
  arrange(desc(post_mean))

# Compare naive sample rates vs EB posterior means
compare_df <- team_stats %>%
  mutate(sample_rate = k/n) %>%
  left_join(team_post_eb %>% select(team, post_mean), by = "team") %>%
  pivot_longer(cols = c(sample_rate, post_mean), names_to = "type", values_to = "p")

ggplot(compare_df, aes(p, fill = type)) +
  geom_histogram(bins = 20, alpha = 0.6, position = "identity") +
  labs(
    title = "Shrinkage Effect: Sample Rates vs EB Posterior Means",
    x = "Probability",
    y = "Count"
  )

Shrinkage reduces overreaction to small samples—one of the most common reasons betting systems look great in backtests and fail live.

8. Bayesian Logistic Models: Match Win Probabilities from Features

Team strength depends on context: home advantage, rest, injuries, travel, matchups, and more. Logistic regression is a natural model for win probabilities. Bayesian logistic regression adds regularization and uncertainty.

We’ll create a few simple features from our toy data. In real projects you’d engineer better predictors.

# Simple feature engineering (toy)
# Use rolling win rate as a proxy strength feature (in real data: Elo, xG, etc.)
matches_feat <- matches %>%
  arrange(date) %>%
  group_by(team) %>%
  mutate(
    games_played = row_number() - 1,
    roll_winrate = ifelse(games_played == 0, NA_real_,
                          cummean(lag(result)))
  ) %>%
  ungroup() %>%
  mutate(
    roll_winrate = ifelse(is.na(roll_winrate), mean(result), roll_winrate),
    # center features
    roll_winrate_c = roll_winrate - mean(roll_winrate),
    home_c = home - mean(home)
  )

# Train-test split by time (avoid leakage)
cut_date <- quantile(matches_feat$date, 0.8)
train <- matches_feat %>% filter(date <= cut_date)
test  <- matches_feat %>% filter(date >  cut_date)

nrow(train); nrow(test)
# Bayesian logistic regression
# Note: rstanarm can be slower. For a big dataset, adjust iter/chains or use variational inference elsewhere.
fit_bayes_logit <- stan_glm(
  result ~ home_c + roll_winrate_c,
  data = train,
  family = binomial(link = "logit"),
  prior = normal(location = 0, scale = 1, autoscale = TRUE),
  prior_intercept = normal(0, 2.5),
  chains = 2,
  iter = 800,
  refresh = 0
)

print(fit_bayes_logit)
# Posterior predictive probabilities on test
p_test <- posterior_linpred(fit_bayes_logit, newdata = test, transform = TRUE)
# p_test is draws x observations; take posterior mean probability
p_hat <- colMeans(p_test)

test_pred <- test %>%
  mutate(p_model = p_hat)

test_pred %>%
  summarise(
    mean_p = mean(p_model),
    mean_y = mean(result)
  )

9. Ratings as Priors: Elo with Uncertainty and Bayesian Flavor

Elo is a practical rating system that translates team strength differences into win probabilities. It’s not “Bayesian” in a strict sense, but it behaves like a sequential updating scheme and can be combined with Bayesian ideas.

9.1 Simple Elo implementation in R

elo_update <- function(Ra, Rb, y, k = 20) {
  # y: 1 if A wins, 0 if B wins
  Ea <- 1 / (1 + 10^((Rb - Ra)/400))
  Ra_new <- Ra + k*(y - Ea)
  Rb_new <- Rb + k*((1 - y) - (1 - Ea))
  c(Ra = Ra_new, Rb = Rb_new, Ea = Ea)
}

run_elo <- function(df, init = 1500, k = 20) {
  teams <- sort(unique(c(df$team, df$opponent)))
  R <- setNames(rep(init, length(teams)), teams)

  out <- vector("list", nrow(df))
  for (i in seq_len(nrow(df))) {
    a <- df$team[i]
    b <- df$opponent[i]
    y <- df$result[i]

    upd <- elo_update(R[a], R[b], y, k = k)
    R[a] <- upd["Ra"]
    R[b] <- upd["Rb"]

    out[[i]] <- tibble(
      date = df$date[i],
      team = a,
      opponent = b,
      result = y,
      elo_team = R[a],
      elo_opp = R[b],
      p_elo = upd["Ea"]
    )
  }
  bind_rows(out)
}

elo_df <- run_elo(matches %>% arrange(date), k = 18)
elo_df %>% slice(1:5)
# Join Elo features back to matches
matches_elo <- matches %>%
  arrange(date) %>%
  mutate(row_id = row_number()) %>%
  left_join(
    elo_df %>% mutate(row_id = row_number()) %>% select(row_id, p_elo),
    by = "row_id"
  )

summary(matches_elo$p_elo)

Elo probabilities can serve as baseline priors or features inside a Bayesian model. The key is to evaluate calibration and stability.

10. Score Models: Poisson and Skellam for Totals and Spreads

If you bet totals (Over/Under) or spreads, it’s helpful to model the score distribution. A common starting point is Poisson scoring for each team: HomeGoals ~ Poisson(lambda_home), AwayGoals ~ Poisson(lambda_away). The goal difference follows a Skellam distribution (difference of two Poissons).

We’ll simulate goals here to demonstrate the workflow. Replace this with real score data.

# Add synthetic goals to the toy dataset (for demonstration)
add_goals <- function(df) {
  # tie lambda to latent p_true to create plausible relationship
  # this is just for tutorial purposes
  base <- 1.2
  spread <- (df$p_true - 0.5) * 2.0

  lambda_home <- pmax(0.2, base + 0.4 + spread)
  lambda_away <- pmax(0.2, base - 0.1 - spread)

  df %>%
    mutate(
      goals_team = rpois(n(), lambda_home),
      goals_opp  = rpois(n(), lambda_away),
      total_goals = goals_team + goals_opp,
      goal_diff = goals_team - goals_opp
    )
}

matches_goals <- add_goals(matches)
summary(matches_goals$total_goals)

10.1 Fit simple Poisson regression for scoring rates

# Poisson model for "team goals" using simple predictors
# In real data: team/opponent effects, home advantage, pace, etc.
train_g <- matches_goals %>% filter(date <= cut_date)
test_g  <- matches_goals %>% filter(date >  cut_date)

fit_pois <- glm(
  goals_team ~ home + p_true,
  data = train_g,
  family = poisson()
)

summary(fit_pois)
# Predict lambdas on test
lambda_hat <- predict(fit_pois, newdata = test_g, type = "response")

test_g2 <- test_g %>% mutate(lambda_team = lambda_hat)

summary(test_g2$lambda_team)

10.2 Convert score distribution into market probabilities (totals/spreads)

To price a total (e.g., Over 2.5), you need P(Total > 2.5). If we model both sides as Poisson and independent, total is Poisson with rate lambda_total = lambda_team + lambda_opp. Here we only have a team-goals model; we’ll create a companion lambda for opponent for the demo.

# Companion opponent lambda (toy)
fit_pois_opp <- glm(
  goals_opp ~ home + p_true,
  data = train_g,
  family = poisson()
)

test_g2 <- test_g2 %>%
  mutate(lambda_opp = predict(fit_pois_opp, newdata = test_g2, type = "response"),
         lambda_total = lambda_team + lambda_opp)

# Price Over 2.5
p_over_25 <- 1 - ppois(2, lambda = test_g2$lambda_total)
summary(p_over_25)
# Example: bet Over 2.5 at decimal odds 1.95
d_over <- 1.95
ev_over <- ev_decimal(p_over_25, d_over)
summary(ev_over)
mean(ev_over > 0)

11. Probability Calibration: Brier, Log Loss, and Reliability Curves

A betting system needs probabilities that mean what they say. If you claim 60% and you win 60% over time, that’s calibration. Calibration is one of the best predictors that your edge is real rather than backtest noise.

11.1 Brier score and log loss

brier_score <- function(y, p) mean((y - p)^2)

log_loss <- function(y, p, eps = 1e-12) {
  p <- pmin(pmax(p, eps), 1 - eps)
  -mean(y*log(p) + (1-y)*log(1-p))
}

brier <- brier_score(test_pred$result, test_pred$p_model)
ll <- log_loss(test_pred$result, test_pred$p_model)
c(brier = brier, logloss = ll)

11.2 Reliability curve (calibration plot)

reliability_curve <- function(y, p, bins = 10) {
  tibble(y = y, p = p) %>%
    mutate(bin = ntile(p, bins)) %>%
    group_by(bin) %>%
    summarise(
      p_mean = mean(p),
      y_mean = mean(y),
      n = n(),
      .groups = "drop"
    )
}

rc <- reliability_curve(test_pred$result, test_pred$p_model, bins = 12)

ggplot(rc, aes(p_mean, y_mean, size = n)) +
  geom_point(alpha = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  labs(
    title = "Reliability Curve (Calibration)",
    x = "Predicted probability (bin mean)",
    y = "Observed frequency",
    size = "Count"
  )

If your points drift above the diagonal, you’re underconfident; below it, overconfident. Both can destroy Kelly sizing.

12. Kelly Logic: Optimal Growth, Fractional Kelly, and Practical Constraints

Kelly sizing maximizes long-run logarithmic bankroll growth if your probabilities are correct. In real life, probabilities are noisy—so full Kelly can be too aggressive. Fractional Kelly is often the sweet spot.

12.1 Kelly fraction for decimal odds

p <- 0.55
d <- 2.05

f_k <- kelly_fraction_decimal(p, d)
f_k
# Fractional Kelly (e.g., half Kelly)
fractional_kelly <- function(f, frac = 0.5) frac * f

f_half <- fractional_kelly(f_k, 0.5)

# Apply caps (never bet negative Kelly; cap max stake)
f_bet <- clip(f_half, 0, 0.05)
f_bet

12.2 Kelly with uncertainty (posterior integration)

A Bayesian twist: rather than plug in a point estimate p, integrate over uncertainty by sampling from the posterior. This yields a distribution for Kelly fraction and lets you bet conservatively (e.g., on a lower quantile).

kelly_draws_decimal <- function(p_draws, d) {
  b <- d - 1
  q <- 1 - p_draws
  (b*p_draws - q)/b
}

# Example using the earlier Beta posterior (pick some a,b)
a <- 25; b <- 20
p_draws <- rbeta(50000, a, b)
d <- 1.95

f_draws <- kelly_draws_decimal(p_draws, d)

quantile(f_draws, c(0.05, 0.25, 0.5, 0.75, 0.95))
# Conservative staking: use 25th percentile of Kelly, then apply fraction + cap
f_cons <- quantile(f_draws, 0.25)
f_bet <- clip(0.5 * f_cons, 0, 0.03)
f_bet

This is one of the most practical Bayesian risk controls: uncertainty shrinks position sizes automatically.

13. Portfolio Betting: Correlation, Market Clustering, and Exposure Controls

Betting is rarely a sequence of independent wagers. Same league, same teams, same injury news—your bets can be correlated. Correlation increases drawdowns and makes naive Kelly sizing too aggressive.

13.1 A simple exposure constraint system

# Imagine we have candidate bets with:
# date, market, team, p_model, odds_decimal, kelly_f

make_candidates <- function(df, p_col = "p_model", odds_col = "odds_decimal") {
  df %>%
    transmute(
      date,
      team,
      opponent,
      p = .data[[p_col]],
      d = .data[[odds_col]],
      ev = ev_decimal(p, d),
      f_kelly = kelly_fraction_decimal(p, d)
    ) %>%
    mutate(
      f_half = 0.5 * f_kelly,
      f_bet = clip(f_half, 0, 0.03)
    )
}

# For demo, use test_pred probabilities against the odds in matches
candidates <- test_pred %>%
  mutate(p_model = p_model, odds_decimal = odds_decimal) %>%
  make_candidates()

summary(candidates$f_bet)
summary(candidates$ev)
# Portfolio constraints:
# - Max stake per day
# - Max stake per team
# - Only bet when EV > 0 and Kelly > 0

apply_constraints <- function(df, max_daily = 0.08, max_team = 0.05) {
  df %>%
    filter(ev > 0, f_bet > 0) %>%
    arrange(desc(ev)) %>%
    group_by(date) %>%
    mutate(
      daily_cum = cumsum(f_bet),
      f_bet_day = ifelse(daily_cum <= max_daily, f_bet, 0)
    ) %>%
    ungroup() %>%
    group_by(team) %>%
    mutate(
      team_cum = cumsum(f_bet_day),
      f_bet_final = ifelse(team_cum <= max_team, f_bet_day, 0)
    ) %>%
    ungroup()
}

bets_planned <- apply_constraints(candidates, max_daily = 0.08, max_team = 0.05)

bets_planned %>%
  summarise(
    n_candidates = n(),
    n_bets = sum(f_bet_final > 0),
    avg_stake = mean(f_bet_final[f_bet_final > 0]),
    total_stake = sum(f_bet_final)
  )

This looks simple, but constraints are often the difference between a backtest hero and a live survivor.

14. Concrete Betting Strategies: From Simple Rules to Model-Driven Systems

Below are practical strategies you can implement. None is “magic.” The point is to define rules, measure them, and iterate.

14.1 Strategy A: Simple value betting (probability edge threshold)

strategy_value_threshold <- function(df, edge_min = 0.02) {
  df %>%
    mutate(
      p_be = 1/d,
      edge = p - p_be
    ) %>%
    filter(edge >= edge_min) %>%
    mutate(
      f_bet = clip(0.5 * kelly_fraction_decimal(p, d), 0, 0.03)
    )
}

bets_A <- candidates %>% strategy_value_threshold(edge_min = 0.02)
summary(bets_A$f_bet)

14.2 Strategy B: Bayesian “probability of positive EV” gate

If you can produce a posterior distribution for p, you can compute P(EV > 0). We’ll demonstrate using a Beta posterior for each team’s win rate as a crude uncertainty model.

# Build EB Beta posteriors from the training window
team_train <- matches %>% filter(date <= cut_date) %>%
  group_by(team) %>%
  summarise(n = n(), k = sum(result), .groups = "drop")

p_hat_train <- team_train$k / team_train$n
prior_ab <- fit_beta_mom(p_hat_train, w = team_train$n)
a0 <- prior_ab["a"]; b0 <- prior_ab["b"]

team_params <- team_train %>%
  mutate(a = a0 + k, b = b0 + (n - k)) %>%
  select(team, a, b)

# Attach team posterior params to candidates
cand_B <- candidates %>%
  left_join(team_params, by = "team")

prob_ev_positive_beta <- function(a, b, d, n_draws = 10000) {
  p_draws <- rbeta(n_draws, a, b)
  mean(ev_decimal(p_draws, d) > 0)
}

cand_B2 <- cand_B %>%
  mutate(
    p_ev_pos = pmap_dbl(list(a, b, d), ~prob_ev_positive_beta(..1, ..2, ..3, n_draws = 6000)),
    f_bet = clip(0.5 * kelly_fraction_decimal(p, d), 0, 0.03)
  )

bets_B <- cand_B2 %>% filter(p_ev_pos >= 0.65, ev > 0, f_bet > 0)

bets_B %>%
  summarise(n_bets = n(), mean_p_ev_pos = mean(p_ev_pos), mean_ev = mean(ev), mean_f = mean(f_bet))

14.3 Strategy C: Calibrated model + fractional Kelly + caps

# If you have a calibrated model probability p_model,
# apply fractional Kelly and cap. Also: require minimal EV.
strategy_fractional_kelly <- function(df, ev_min = 0.01, frac = 0.5, cap = 0.03) {
  df %>%
    mutate(
      f_k = kelly_fraction_decimal(p, d),
      f_bet = clip(frac * f_k, 0, cap)
    ) %>%
    filter(ev >= ev_min, f_bet > 0)
}

bets_C <- candidates %>% strategy_fractional_kelly(ev_min = 0.01, frac = 0.35, cap = 0.025)
bets_C %>% summarise(n_bets = n(), mean_ev = mean(ev), mean_f = mean(f_bet))

15. Backtesting in R: Walk-forward, Leakage Prevention, and Metrics

Most betting “edges” vanish when you backtest properly. Your backtest must be time-respecting (walk-forward), must avoid target leakage, and must include realistic constraints (limits, stake caps, minimum odds, etc.).

15.1 A minimal bet settlement engine

settle_bets <- function(df, bankroll0 = 1.0) {
  # df must have: result (1/0), d (decimal odds), f_bet (fraction of bankroll)
  bankroll <- bankroll0
  out <- vector("list", nrow(df))

  for (i in seq_len(nrow(df))) {
    stake <- bankroll * df$f_bet[i]
    if (stake <= 0) {
      out[[i]] <- tibble(i = i, bankroll = bankroll, stake = 0, pnl = 0)
      next
    }

    pnl <- if (df$result[i] == 1) stake * (df$d[i] - 1) else -stake
    bankroll <- bankroll + pnl

    out[[i]] <- tibble(i = i, bankroll = bankroll, stake = stake, pnl = pnl)
  }

  bind_cols(df, bind_rows(out))
}

# Example: settle Strategy C bets (in arbitrary order; ideally sort by date)
bt_C <- bets_C %>%
  mutate(result = test_pred$result[match(paste(date, team, opponent), paste(test_pred$date, test_pred$team, test_pred$opponent))]) %>%
  arrange(date) %>%
  select(date, team, opponent, result, d, f_bet)

bt_C2 <- settle_bets(bt_C, bankroll0 = 1.0)

tail(bt_C2)

15.2 Performance metrics

compute_metrics <- function(df) {
  df %>%
    summarise(
      n_bets = sum(stake > 0),
      final_bankroll = last(bankroll),
      total_pnl = sum(pnl),
      roi_on_staked = ifelse(sum(stake) > 0, sum(pnl)/sum(stake), NA_real_),
      hit_rate = mean(result[stake > 0] == 1),
      avg_odds = mean(d[stake > 0]),
      max_drawdown = {
        peak <- cummax(bankroll)
        dd <- bankroll/peak - 1
        min(dd)
      }
    )
}

compute_metrics(bt_C2)
# Plot bankroll curve
ggplot(bt_C2, aes(date, bankroll)) +
  geom_line() +
  labs(
    title = "Backtest Bankroll Curve (Strategy C)",
    x = "Date",
    y = "Bankroll"
  )

15.3 Walk-forward template

Below is a simple walk-forward skeleton: refit model on an expanding window, predict the next chunk, generate bets, settle, repeat. Replace the “model” section with your real pipeline.

walk_forward <- function(df, initial_frac = 0.6, step_frac = 0.1) {
  df <- df %>% arrange(date)
  n <- nrow(df)

  idx0 <- floor(initial_frac * n)
  step <- floor(step_frac * n)

  start <- idx0
  bankroll <- 1.0
  all_bets <- list()

  while (start < n) {
    train_idx <- 1:start
    test_idx <- (start + 1):min(n, start + step)

    train <- df[train_idx, ]
    test  <- df[test_idx, ]

    # --- Model (placeholder) ---
    # simple baseline: probability = historical mean adjusted by home
    p0 <- mean(train$result)
    p_hat <- clip(p0 + 0.03*(test$home - mean(train$home)), 0.02, 0.98)

    cand <- test %>%
      transmute(
        date, team, opponent,
        result,
        p = p_hat,
        d = odds_decimal
      ) %>%
      mutate(
        ev = ev_decimal(p, d),
        f_bet = clip(0.35 * kelly_fraction_decimal(p, d), 0, 0.02)
      ) %>%
      filter(ev > 0.01, f_bet > 0)

    # Settle on this chunk
    if (nrow(cand) > 0) {
      settled <- settle_bets(cand, bankroll0 = bankroll)
      bankroll <- last(settled$bankroll)
      all_bets[[length(all_bets) + 1]] <- settled
    }

    start <- start + step
  }

  bind_rows(all_bets)
}

wf_bt <- walk_forward(matches)
compute_metrics(wf_bt)

16. Monte Carlo Stress Tests: Drawdowns, Ruin, and Regret

Backtests are one path. Stress tests are another. Monte Carlo simulation answers: “If my estimated edge is real but noisy, what drawdowns should I expect?”

16.1 Simulate a season of bets with uncertainty around p

simulate_betting_path <- function(n_bets = 500, p = 0.54, d = 1.95, f = 0.01, bankroll0 = 1) {
  bankroll <- bankroll0
  for (i in 1:n_bets) {
    stake <- bankroll * f
    win <- rbinom(1, 1, p)
    pnl <- if (win == 1) stake*(d - 1) else -stake
    bankroll <- bankroll + pnl
  }
  bankroll
}

# Many simulations
sim_paths <- replicate(5000, simulate_betting_path(n_bets = 400, p = 0.53, d = 1.95, f = 0.012))
quantile(sim_paths, c(0.05, 0.25, 0.5, 0.75, 0.95))

16.2 Incorporate probability error (model risk)

Your p_hat is not the true p. Model risk is often larger than market risk. We simulate true p drifting around your estimate.

simulate_with_model_error <- function(n_bets = 500, p_hat = 0.54, p_sd = 0.03, d = 1.95, f = 0.01, bankroll0 = 1) {
  bankroll <- bankroll0
  for (i in 1:n_bets) {
    p_true <- clip(rnorm(1, p_hat, p_sd), 0.01, 0.99)
    stake <- bankroll * f
    win <- rbinom(1, 1, p_true)
    pnl <- if (win == 1) stake*(d - 1) else -stake
    bankroll <- bankroll + pnl
  }
  bankroll
}

sim2 <- replicate(5000, simulate_with_model_error(n_bets = 400, p_hat = 0.53, p_sd = 0.05, d = 1.95, f = 0.012))
quantile(sim2, c(0.05, 0.25, 0.5, 0.75, 0.95))
# Compare distributions
df_sim <- tibble(
  bankroll = c(sim_paths, sim2),
  scenario = rep(c("no_model_error", "with_model_error"), each = length(sim_paths))
)

ggplot(df_sim, aes(bankroll)) +
  geom_histogram(bins = 50) +
  facet_wrap(~scenario, scales = "free_y") +
  labs(
    title = "Monte Carlo Bankroll Distribution",
    x = "Final bankroll",
    y = "Count"
  )

If model error collapses your distribution, reduce Kelly fraction, add EV thresholds, and improve calibration.

17. A Go-Live Checklist: Operational and Statistical Guardrails

Time split always: train only on past, predict future.
Log everything: inputs, odds, timestamps, model version, decisions, stakes.
Calibrate probabilities: reliability curves, Brier, log loss.
Respect uncertainty: use posterior bands or conservative probability adjustments.
Size bets defensively: fractional Kelly + caps + daily exposure limits.
Account for correlation: cluster exposures by league/team/news.
Stress-test: Monte Carlo with model error, not just “true p”.
Beware selection bias: strategies discovered by brute force often overfit.
Prefer stability: smaller edge + stable behavior beats fragile high edge.

Table of Contents