使用SymPy求解多方程多变量计算三方真实赔率的问题排查

Let's walk through what's going wrong with your current code and how to adjust your approach to get the true probabilities you're expecting.

First: Understand the Core Problem #

You're trying to calculate the "true" probabilities behind bookmaker odds (removing the bookmaker's margin, or "vig"). Your initial approach uses a power transformation, but your system of equations has redundant and conflicting constraints that are causing SymPy to return a parameterized, unhelpful solution.

Why Your Current Code Fails #

Let's break down your equations:

• You define x = 1/odds[0], y = 1/odds[1], z = 1/odds[2] — these are the "implied probabilities" from the odds, which sum to more than 1 (since bookmakers add a margin to guarantee profit).
• Then you set p^n = x, q^n = y, r^n = z — which means p = x^(1/n), q = y^(1/n), r = z^(1/n) by definition.
• Adding p + q + r = 1 is the key constraint for true probabilities, but your other equations (p^n + q^n + r^n = x + y + z and x^(1/n) + y^(1/n) + z^(1/n) = p + q + r) are redundant: they just restate definitions or known values.

SymPy can't solve this because you're essentially asking it to solve for n with only one meaningful constraint (x^(1/n) + y^(1/n) + z^(1/n) = 1), but your other equations don't add new information. The parameterized result {(-q - r + 1, q, r, n)} is SymPy's way of saying "p is determined by q and r via the last equation, but q, r, n can be any values that fit the other redundant rules".

Two Valid Approaches to Calculate True Probabilities #

Approach 1: The Standard "Vig Removal" Method (Simplest) #

The most common way to get true probabilities is to scale the implied probabilities by the payout ratio (the inverse of the bookmaker's margin):

1. Calculate the sum of implied probabilities: total_implied = sum(1/odd for odd in odds)
2. The payout ratio is k = 1 / total_implied (this is how much of each bet the bookmaker actually pays out)
3. True probabilities are p = (1/odds[0]) * k, q = (1/odds[1]) * k, r = (1/odds[2]) * k

This ensures p + q + r = 1, and each true probability is slightly lower than the implied probability (which matches your expectation of "values slightly below x, y, z").

Here's the SymPy code for this:

from sympy import Rational

odds = [1.41, 5.09, 8.4]
# Convert odds to rational implied probabilities for better precision
x = Rational(1, odds[0])
y = Rational(1, odds[1])
z = Rational(1, odds[2])

total_implied = x + y + z
k = 1 / total_implied

true_p = x * k
true_q = y * k
true_r = z * k

print(f"True probabilities: {true_p.evalf()}, {true_q.evalf()}, {true_r.evalf()}")
print(f"Sum: {(true_p + true_q + true_r).evalf()}") # Will equal 1

Approach 2: Power Transformation (Your Original Idea) #

If you specifically want to find a power n such that p = x^(1/n), q = y^(1/n), r = z^(1/n) and p + q + r = 1, you need to simplify your system to focus on solving for n, then compute p, q, r.

Since this equation has no analytical solution, we'll use SymPy's numerical solver nsolve:

from sympy import symbols, Eq, nsolve

n = symbols('n', real=True, positive=True)
odds = [1.41, 5.09, 8.4]
x = 1/odds[0]
y = 1/odds[1]
z = 1/odds[2]

# The only meaningful constraint: sum of transformed probabilities equals 1
equation = Eq(x**(1/n) + y**(1/n) + z**(1/n), 1)

# Solve for n numerically (start with an initial guess, e.g., 1.2)
n_solution = nsolve(equation, n, 1.2)

# Calculate true probabilities
p = x**(1/n_solution)
q = y**(1/n_solution)
r = z**(1/n_solution)

print(f"n = {n_solution.evalf()}")
print(f"True probabilities: {p.evalf()}, {q.evalf()}, {r.evalf()}")
print(f"Sum: {(p + q + r).evalf()}") # Will equal 1

Adding Inequality Constraints (If Needed) #

If you want to enforce x > p, y > q, z > r (which holds in both approaches since the bookmaker's margin makes implied probabilities higher than true probabilities):

• For Approach 1: Since k = 1 / total_implied and total_implied > 1 (due to vig), k < 1, so true_p = x*k < x automatically. No need for extra constraints.
• For Approach 2: Since n > 1 (because x^(1/n) < x when n>1 and 0<x<1), you can restrict n to positive values greater than 1 in SymPy by adding positive=True and using a solver that respects bounds, or verify after solving that the inequalities hold.

Final Notes #

• The standard vig removal method is the industry standard because it's straightforward and mathematically sound for adjusting bookmaker odds to true probabilities.
• Your power transformation approach is a valid alternative but is less common; it assumes a specific functional form for scaling down the implied probabilities.

内容的提问来源于stack exchange，提问作者Max Mustermann