构建指定支持域的截断Beta分布：代码报错与均值偏差求助

问题需求 #

需要指定一个支持域为[0.1, 0.4]的截断Beta分布，要求下界处的概率密度远高于0附近。

第一次尝试与报错 #

代码实现 #

import numpy as np
import pandas as pd
import plotly.express as px
from scipy import stats

def get_ab(mean, stdev):
    kappa = (mean*(1-mean)/stdev**2 - 1)
    return mean*kappa, (1-mean)*kappa

class truncated_beta(stats.rv_continuous):
    def _pdf(self, x, alpha, beta, a, b):
        return stats.beta.pdf(x, alpha, beta) / (stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))
    def _cdf(self, x, alpha, beta, a, b):
        return (stats.beta.cdf(x, alpha, beta) - stats.beta.cdf(a, alpha, beta)) / (stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))
   
def plot_beta_distr(mu, stdev, lb, ub):
    alpha_, beta_ = get_ab((mu-lb)/ub, stdev/ub)

    dist = truncated_beta(a=lb, b=ub, alpha = alpha_, beta = beta_, name='truncated_beta')

    x = np.linspace(lb, ub, 100000)
    y = dist.pdf(x, alpha = alpha_, beta = beta_, a = lb, b = ub)

    fig = px.line(x = x, y = y)
    return fig

mu = 0.02
stdev = 0.005
lb = 0.1
ub = 0.4

plot_beta_distr(mu, stdev, lb, ub)

报错信息 #

TypeError: __init__() got an unexpected keyword argument 'alpha'

问题排查 #

报错有两个核心原因：

1. 指定的mu=0.02超出了截断区间[0.1, 0.4]的范围，导致后续参数计算逻辑完全失效
2. 自定义truncated_beta类时，错误地将alpha、beta作为实例化参数传入，而stats.rv_continuous的__init__方法并不接受这类参数

第二次尝试与问题 #

代码实现 #

mean = .35
std = .1

lb = 0.1
ub = 0.5

def get_ab(mu, stdev):
    kappa = (mu*(1-mu)/stdev**2 - 1)
    return mu*kappa, (1-mu)*kappa

alpha_, beta_ = get_ab((mean - lb) / ub, std/ub)

norm = stats.beta.cdf(ub, alpha_, beta_) - stats.beta.cdf(lb, alpha_, beta_)

yr=stats.uniform.rvs(size=1000000)*norm+stats.beta.cdf(lb, alpha_, beta_)
xr=stats.beta.ppf(yr, alpha_, beta_)

xr.min()  # 0.100
xr.max()  # 0.4999
xr.std()  # 0.104
xr.mean() # 0.341

存在问题 #

样本的最小值、最大值、标准差符合预期，但均值与设定的.35不符。核心错误是原始Beta分布的参数转换逻辑错误：直接对均值做(mean - lb)/ub的线性缩放并不适用于截断Beta分布，截断后的分布均值无法通过简单缩放原始Beta分布均值得到。

正确实现方案 #

思路说明 #

要生成支持域为[a, b]的截断Beta分布，需明确两点：

1. 原始Beta分布定义在[0,1]区间，截断到[a,b]后，概率密度是原始Beta密度在[a,b]上的归一化结果
2. 不能通过简单缩放原始Beta分布的均值来匹配截断后的目标均值，必须使用矩匹配法结合数值优化求解合适的alpha和beta参数

代码实现 #

import numpy as np
from scipy import stats
from scipy.optimize import root_scalar, root

# 目标参数
target_mean = 0.35
target_std = 0.1
lb = 0.1
ub = 0.4

def truncated_beta_mean(alpha, beta, a, b):
    """计算[0,1]上Beta分布截断到[a,b]后的均值"""
    mu0 = alpha / (alpha + beta)
    cdf_a = stats.beta.cdf(a, alpha, beta)
    cdf_b = stats.beta.cdf(b, alpha, beta)
    pdf_a = stats.beta.pdf(a, alpha, beta)
    pdf_b = stats.beta.pdf(b, alpha, beta)
    numerator = mu0*(cdf_b - cdf_a) - (b*pdf_b - a*pdf_a)/(alpha+beta)
    return numerator / (cdf_b - cdf_a)

def truncated_beta_var(alpha, beta, a, b):
    """计算[0,1]上Beta分布截断到[a,b]后的方差"""
    mu_trunc = truncated_beta_mean(alpha, beta, a, b)
    second_moment = ( (alpha*(alpha+1))/((alpha+beta)*(alpha+beta+1)) )*(stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))
    second_moment -= (b**2*stats.beta.pdf(b, alpha, beta) - a**2*stats.beta.pdf(a, alpha, beta))/(alpha+beta)
    second_moment /= (stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))
    return second_moment - mu_trunc**2

def solve_alpha_beta(target_mean, target_std, a, b):
    """通过数值优化求解符合目标均值和方差的alpha、beta"""
    target_var = target_std**2
    
    # 先估计初始alpha值
    def mean_error(alpha):
        return truncated_beta_mean(alpha, 1, a, b) - target_mean
    res = root_scalar(mean_error, bracket=[0.1, 100], method='brentq')
    alpha_init = res.root
    
    # 再估计初始beta值
    def var_error(beta):
        return truncated_beta_var(alpha_init, beta, a, b) - target_var
    res = root_scalar(var_error, bracket=[0.1, 100], method='brentq')
    beta_init = res.root
    
    # 联合优化alpha和beta，同时匹配均值和方差
    def equations(params):
        alpha, beta = params
        mu = truncated_beta_mean(alpha, beta, a, b)
        var = truncated_beta_var(alpha, beta, a, b)
        return [mu - target_mean, var - target_var]
    
    opt_res = root(equations, [alpha_init, beta_init])
    return opt_res.x[0], opt_res.x[1]

# 求解最优参数
alpha_opt, beta_opt = solve_alpha_beta(target_mean, target_std, lb, ub)
print(f"最优alpha: {alpha_opt:.4f}, 最优beta: {beta_opt:.4f}")

# 生成截断Beta分布样本
cdf_lb = stats.beta.cdf(lb, alpha_opt, beta_opt)
cdf_ub = stats.beta.cdf(ub, alpha_opt, beta_opt)
yr = stats.uniform.rvs(size=1000000)*(cdf_ub - cdf_lb) + cdf_lb
xr = stats.beta.ppf(yr, alpha_opt, beta_opt)

# 验证样本统计量
print(f"样本均值: {xr.mean():.4f}, 目标均值: {target_mean}")
print(f"样本标准差: {xr.std():.4f}, 目标标准差: {target_std}")
print(f"样本最小值: {xr.min():.4f}, 下界: {lb}")
print(f"样本最大值: {xr.max():.4f}, 上界: {ub}")

关键说明 #

1. 实现了截断Beta分布的均值和方差精确计算公式，这是矩匹配的核心依据
2. 通过数值优化方法迭代求解alpha和beta，确保截断后的分布完全匹配目标均值和方差
3. 样本生成逻辑保持不变，但使用优化后的参数，最终样本统计量会与目标值一致

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