使用SymPy求解Sortino比率半标准差积分解析解报错求助

问题背景 #

需计算Sortino比率对应的半标准差积分，积分表达式为：
$$\sqrt{\int_{-\infty}^{B} (B - x)^2 f(x) dx}$$
（注：积分区间内$B-x\geq0$，$\sqrt{(B-x)^2}$可简化为$B-x$）

原代码与报错 #

原代码 #

import sympy as sy 
c = sy.Symbol('B')
x = sy.Symbol('x') 
a= sy.Symbol('-sy.oo')

x = sy.Symbol('x')
f= sy.Symbol('f')
d= sy.Symbol('d')

f1 = sy.integrate(sy.sqrt((c-x)**2* f(x)*dx),(x, -sy.oo, c))

报错信息 #

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/hc/13tll2g535x553nmqlfcyztr0000gn/T/ipykernel_70280/1712613650.py in <module>
----> 1 f1 = sy.integrate(sy.sqrt((c-x)**2* f(x)*dx),(x, -sy.oo, 0.2))

TypeError: 'Symbol' object is not callable

报错原因与修正方案 #

核心错误点 #

1. Symbol被当作函数调用：f被定义为SymPy符号，却以f(x)形式调用，触发"Symbol不可调用"错误。
2. 冗余的dx符号：SymPy的integrate函数无需手动传入dx，积分变量由第二个参数元组指定。
3. 错误定义负无穷：-sy.oo是SymPy内置的负无穷常量，无需定义为Symbol。

修正后的代码（以正态分布为例） #

Sortino比率通常假设收益服从正态分布，以下是针对该分布的解析解计算代码：

import sympy as sy

# 定义符号变量
x = sy.Symbol('x')
B = sy.Symbol('B')  # 阈值收益率
mu = sy.Symbol('mu')  # 收益均值
sigma = sy.Symbol('sigma')  # 收益标准差

# 定义正态分布概率密度函数f(x)
f_x = (1/(sigma * sy.sqrt(2 * sy.pi))) * sy.exp(-(x - mu)**2/(2 * sigma**2))

# 简化被积函数，计算半标准差（积分结果的平方根）
integrand = (B - x)**2 * f_x
sortino_semistd = sy.sqrt(sy.integrate(integrand, (x, -sy.oo, B)))

# 打印化简后的解析解
print(sy.simplify(sortino_semistd))

解析解说明 #

• 仅当f(x)为特定已知分布（如正态、对数正态）时，SymPy才能生成解析解。若f(x)是无具体形式的通用概率密度函数，无法得到通用解析解。
• 正态分布下的半标准差解析解化简后为：
  $$\sigma \sqrt{\frac{1}{2} \left(1 + \text{erf}\left(\frac{B - \mu}{\sigma\sqrt{2}}\right)\right)} + \frac{(\mu - B)}{\sqrt{2\pi}} e^{-\frac{(B - \mu)^2}{2\sigma2}}$$
  （最终半标准差为上述表达式的平方根）

数值积分备选方案 #

若无法获取解析解（如自定义分布），可使用数值积分计算：

from scipy.integrate import quad
import numpy as np

# 代入具体参数值
B_val = 0.2
mu_val = 0.1
sigma_val = 0.15

# 定义数值被积函数
def integrand_numeric(x):
    pdf = (1/(sigma_val * np.sqrt(2 * np.pi))) * np.exp(-(x - mu_val)**2/(2 * sigma_val**2))
    return (B_val - x)**2 * pdf

# 计算数值积分并求平方根
integral_result, _ = quad(integrand_numeric, -np.inf, B_val)
sortino_semistd_numeric = np.sqrt(integral_result)
print(sortino_semistd_numeric)

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