来瓶青柑普洱吧！

一个基于另一种假设（并且十分概率论）的分析：czm233 的概率论大学习

#Prototype

东方树叶搞活动了，形如：

购买本产品，扫描瓶盖内侧二维码，即有机会赢取：

• $666$ 元红包（中奖率 $0.00008\%$）
• $66$ 元红包（中奖率 $0.005\%$）
• $2$ 元红包（中奖率 $0.5\%$）
• $1$ 元红包（中奖率 $1.0\%$）
• $0.5$ 元红包（中奖率 $16.5\%$）
• 壹元换购（中奖率 $25\%$）

按“能换购则换购”的策略，期望单价是多少？

#Problems

现在考虑一个简化版的问题：

某饮料售价为 $c$，有 $a$ 的概率获得 $r$ 元，有 $b$ 的概率可以花 $w$ 元复购。将进行一次购买和若干次 $w$ 元复购记为一轮，现在研究以下问题：

1. 若希望购买 $N$ 瓶，求期望单价

2. 若希望购买 $M$ 轮，求期望单价

#Solution of 1

既然总要买满 $N$ 瓶，那么复购的唯一意义就是降低花费，也就是等价于若购买下一瓶，则获得 $(c-w)$ 元的补贴。故有：

$$\begin{aligned} E_{1,N}(cost) &= \dfrac{(N - 1)(c - ar - b(c - w)) + (c - ar)}{N} = c-ar - \dfrac{b(N-1)(c-w)}{N}\\ E_{1,\infty}(cost) &= c - ar - b(c-w) \end{aligned}$$

#Solution of 2

先考虑一轮购买。假设购得 $1 + k$ 瓶，记此事件为 $A_{k}$：

$$\begin{aligned} Pr(A_k) &= b^k(1-b)\\ E(cost_k) &= c + kw - \dfrac{ar}{1-b} \end{aligned}$$

为了研究方便，为 $Pr(A_k)$ 构造一个生成函数：

$$f(x) = (1-b)\sum_{k=0}^{\infty} (xb)^k = \dfrac{1-b}{1-xb}\\$$

现在假设有 $M$ 轮购买，共购得 $M + p$ 瓶，记此事件为 $B_{p}$。先研究概率，由 $Pr(A_k)$ 构造生成函数可得：

$$\begin{aligned} Pr(B_p) &= [x^p] \left(\dfrac{1-b}{1-xb}\right)^M\\ E(Bcost_p) &= Mc + pw - \dfrac{Mar}{1-b}=Mu+pw(u=c-\dfrac{ar}{1-b})\\ E_{2,M}(cost) &= \sum_{p = 0}^{\infty} Pr(B_p)\dfrac{E(Bcost_p)}{M + p}\\ &= w + M(u-w)\sum_{p = 0}^{\infty} \dfrac{Pr(B_p)}{M + p} \end{aligned}$$

然后等于我不会了。

#Solution of 2 - rework

首先换一个形式：单价 $c$，复购概率 $p$，复购价格 $b$，共 $k$ 轮。

单价是拿红包价值算出来的，总之有一些方法可以转换回去。

设一共买了 $n$ 瓶，则有：

$$\begin{aligned} E &= \left(\frac{1-p}{p}\right)^k \sum_{n = k}^{\infty} \left(\binom{n-1}{k-1} p^{n} \frac{ck+b(n-k)}{n}\right)\\ &= b + (c-b)k \left( \frac{1-p}{p} \right)^k \left[ \int_0^{\frac{p}{1-p}} \frac{u^{k-1}}{1+u} du \right]\\ &= b + (c-b)k \left( \frac{1-p}{p} \right)^k \left[ (-1)^k \ln(1-p) + \sum_{j=1}^{k-1} (-1)^{k-1-j} \frac{1}{j} \left( \frac{p}{1-p} \right)^j \right]\\ \end{aligned}$$

求一下 $k \to \infty$ 时的收敛值：

$$E = b + (c-b)k \left( \frac{1-p}{p} \right)^k \left[ \int_0^{\frac{p}{1-p}} \frac{u^{k-1}}{1+u} du \right]$$

令 $x = \frac{p}{1-p}$。则 $0 < p < 1$ 意味着 $x > 0$。

我们关注积分 $I_k = \int_0^{x} \frac{u^{k-1}}{1+u} du$ 在 $k \to \infty$ 时的行为。对于这类积分，可以使用拉普拉斯方法来估计其渐近行为。

对于 $x > 0$，当 $k \to \infty$ 时，积分的贡献主要来自于 $u$ 接近上限 $x$ 的区域。通过分部积分法，可以得到渐近展开式：

$$I_k = \int_0^{x} \frac{u^{k-1}}{1+u} du = \frac{x^k}{k(1+x)} + O\left(\frac{x^k}{k^2}\right)$$

这个近似在 $x>0$ 时都成立。将 $I_k$ 的渐近表达式代入 $E$ 中：

$$\begin{aligned} E &= b + (c-b)k \left( \frac{1-p}{p} \right)^k \left[ \frac{1}{k} \left( \frac{p}{1-p} \right)^k \frac{1}{1+\frac{p}{1-p}} + O\left(\frac{1}{k^2} \left(\frac{p}{1-p}\right)^{k+1} \right) \right]\\ &= b + (c-b) \left[ \left( \frac{1-p}{p} \right)^k \left(\frac{p}{1-p}\right)^k \frac{1}{1+\frac{p}{1-p}} + k \cdot O\left(\frac{1}{k^2} \left(\frac{1-p}{p}\right)^k \left(\frac{p}{1-p}\right)^{k+1} \right) \right]\\ &= b + (c-b) \left[ 1-p + k \cdot O\left(\frac{1}{k^2} \frac{p}{1-p} \right) \right] \end{aligned}$$

当 $k \to \infty$ 时，余项 $O\left(\frac{1}{k} \frac{p}{1-p} \right)$ 趋近于 $0$。将简化后的项代回 $E$ 的表达式，并取极限：

$$\lim_{k \to \infty} E = b + (c-b)(1-p)= bp + c(1-p)$$

做完了！问 Gemini 跑个蒙特卡洛试试看，看起来还是挺对的。

import math

def calculate_expected_price(k: int, b: float, c: float, p: float) -> float:
    """
    根据封闭形式的数学公式精确计算k轮购买后的平均单价期望。

    Args:
        k (int): 购买的轮数 (必须为正整数)。
        b (float): 换购价格。
        c (float): 首次购买价格。
        p (float): 获得换购资格的概率 (0 < p < 1)。

    Returns:
        float: 总的平均单价的期望。
    """
    # --- 参数校验 ---
    if not isinstance(k, int) or k < 1:
        raise ValueError("k 必须是正整数")
    if not 0 < p < 1:
        raise ValueError("概率 p 必须在 0 和 1 之间")

    # --- 特殊情况处理 ---
    if b == c:
        return b

    # --- 根据公式计算 ---
    # 第一部分
    result = b

    # 计算公式的第二部分
    log_term = ((-1)**k) * math.log(1 - p)

    sum_term = 0.0
    # 计算求和项 Σ(...)
    # 当 k=1 时, range(1, 1) 为空，此循环不执行，sum_term 为 0
    for j in range(1, k):
        term = ((-1)**(k - 1 - j)) / j * math.pow(p / (1 - p), j)
        sum_term += term

    # 方括号内的总和
    inside_brackets = log_term + sum_term
    
    # 方括号外的系数
    coefficient = (c - b) * k * math.pow((1 - p) / p, k)

    result += coefficient * inside_brackets
    return result

import random

def simulate_average_price(k: int, b: float, c: float, p: float, num_simulations: int = 100000) -> float:
    """
    使用蒙特卡洛方法模拟计算k轮购买后的平均单价期望。

    Args:
        k (int): 购买的轮数 (必须为正整数)。
        b (float): 换购价格。
        c (float): 首次购买价格。
        p (float): 获得换购资格的概率 (0 < p < 1)。
        num_simulations (int): 模拟实验的总次数，次数越多结果越精确。

    Returns:
        float: 总的平均单价的近似期望。
    """
    # --- 参数校验 ---
    if not isinstance(k, int) or k < 1:
        raise ValueError("k 必须是正整数")
    if not 0 < p < 1:
        raise ValueError("概率 p 必须在 0 和 1 之间")

    total_of_average_prices = 0.0

    # 进行大量独立的模拟实验
    for _ in range(num_simulations):
        
        # 在单次实验中，计算k轮购买的总成本和总瓶数
        total_cost_k_rounds = 0
        total_bottles_k_rounds = 0

        for _ in range(k):
            # --- 模拟一轮购买 ---
            # 先以原价c购买一瓶
            round_cost = c
            round_bottles = 1
            
            # 只要抽奖成功(random() < p)，就以换购价b继续购买
            while random.random() < p:
                round_cost += b
                round_bottles += 1
            
            # 累加这一轮的成本和瓶数
            total_cost_k_rounds += round_cost
            total_bottles_k_rounds += round_bottles
            
        # 计算本次实验（k轮购买）的平均单价
        current_average_price = total_cost_k_rounds / total_bottles_k_rounds
        
        # 将本次实验的平均单价累加到总和中
        total_of_average_prices += current_average_price

    # 所有实验的平均单价的平均值，作为期望的近似
    return total_of_average_prices / num_simulations

k = 9
c = 10.0
b = 2.0
p = 0.75

exact_value = calculate_expected_price(k, b, c, p)

simulated_value = simulate_average_price(k, b, c, p, num_simulations=10000000)

print(f"Args: k={k}, c={c}, b={b}, p={p}\n")
print(f"Exact Value: {exact_value:.6f}")
print(f"Monte Carlo: {simulated_value:.6f}")