采用数值二分法计算Hill系数时程序陷入无限循环求助

我正尝试计算两个logistic函数f(x)、g(x)及其复合函数c(x)的Hill系数，并将函数表达式、复合函数、f(x)与g(x)的Hill系数乘积等信息存入DataFrame数据表。计算Hill系数时使用数值二分法求解EC10和EC90，但程序陷入无限循环，强制停止后报错如下：

原代码 #

# Imports
import numpy as np
from numpy import log as ln
from matplotlib import pyplot as plt
from random import randint
import sympy as sym
import pandas as pd
import math


# initializing data
data = {'row': [],
        'f(x)': [],
        'g(x)': [],
        'f(g(x))': [],
        'H_f': [],
        'H_g': [],
        'H_fg': [],
        'Product of H_f and H_g': [],
        'Does this prove hypothesis?': []
        }

df = pd.DataFrame(data)

#True Randomization
def logrand(b):
    ra = randint(1,b)
    ra = float(ra/100)
    ra = 10**ra
    ra = round(ra)
    ra = int(ra)
    return ra

# Two Hill functions
num1 = 10
number = 1
for _ in range(num1):
        # Params
        u = 10
        c1 = logrand(300)
        r1 = randint(1,u)
        k1 = logrand(300)

        c2 = logrand(300)
        r2 = randint(1,u)
        k2 = logrand(300)

        # function layout
        funcf = '{}/({}+e^(-{}x))'.format(c1, k1, r1)
        funcg = '{}/({}+e^(-{}x))'.format(c2, k2, r2)
        funcc = '{}/({}+e^(-{}({}/({}+e^(-{}x)))))'.format(c1, k1, r1, c2, k2, r2)

        # figure layout
        plt.rcParams["figure.figsize"] = [7.50, 3.50]
        plt.rcParams["figure.autolayout"] = True


        # Hill Function for f(x) and g(x) and f(g(x))
        def f(x):
                return c1 / (k1 + np.exp((-1*r1) * x))


        def g(x):
            return c2 / (k2 + np.exp((-1*r2) * x))


        def comp(x):
                return c1 / (k1 + np.exp((-1*r1) * (c2 / (k2 + np.exp((-1*r2) * x)))))



        # EC finder
        def BisectionEC10(fa, a, b):
                c = 1
                x = np.linspace(a, b, 1000)
                ystar = 0.10 * (fa(x).max() - fa(x).min())
                while abs(fa(c) - ystar) > 0.000000001:
                        c = (a + b) / 2
                        if fa(c) - ystar < 0:
                                a = c
                        elif fa(c) - ystar > 0:
                                b = c
                # print('The EC10 of the function is: ',"{0:.15f}".format(c))
                # print('Output of the function when evaluated at the EC10: ',fa(c))
                return c



        def BisectionEC90(fa, a, b):
                c = 1
                x = np.linspace(a, b, 1000)
                ystar = 0.90 * (fa(x).max() - fa(x).min())
                while abs(fa(c) - ystar) > 0.000000001:
                        c = (a + b) / 2
                        if fa(c) - ystar < 0:
                                a = c
                        elif fa(c) - ystar > 0:
                                b = c
                # print('The EC90 of the function is: ',"{0:.15f}".format(c))
                # print('Output of the function when evaluated at the EC90: ',fa(c))
                return c



        # EC90 and EC10 for f(x), g(x) and f(g(x))
        up = 20
        lo = 0
        # x = np.linspace[lo,up,1000]
        x = 1
        # x = sym.symbols('x')

        EC90_1 = BisectionEC90(f, lo, up)
        EC10_1 = BisectionEC10(f, lo, up)

        EC90_2 = BisectionEC90(g, lo, up)
        EC10_2 = BisectionEC10(g, lo, up)

        EC90_3 = BisectionEC90(comp, lo, up)
        EC10_3 = BisectionEC10(comp, lo, up)

        # Hill Coefficient for f(x) and g(x)
        H_1 = ln(81) / (ln(EC90_1 / EC10_1))
        H_1 = round(H_1,4)
        H_2 = ln(81) / (ln(EC90_2 / EC10_2))
        H_2 = round(H_2,4)
        H_3 = ln(81) / (ln(EC90_3 / EC10_3))
        H_3 = round(H_3,4)


        prod = float(H_1.real) * float(H_2.real)
        if prod >= float(H_3.real):
                answer = 'yes'
        else:
                answer = 'no'
        prod = round(prod,4)

        # adding all data to dataframe 2
        data2 = {'row': [number],
                 'f(x)': [funcf],
                 'g(x)': [funcg],
                 'f(g(x))': [funcc],
                 'H_f': [H_1],
                 'H_g': [H_2],
                 'H_fg': [H_3],
                 'Product of H_f and H_g': [prod],
                 'Does this prove hypothesis?': [answer]
                 }
        number = number + 1
        df = df.append(data2, ignore_index=True)

#final dataframe
print(df)
df.to_csv(r'/Users/*****/Desktop/Research/twoarctanfuncs.csv', index = False)

报错信息 #

Traceback (most recent call last):
  File "/Users/*****/Documents/Python/NumericBisectionMethod/venv/twologisticfuncs.py", line 112, in <module>
    EC90_1 = BisectionEC90(f, lo, up)
  File "/Users/*****/Documents/Python/NumericBisectionMethod/venv/twologisticfuncs.py", line 95, in BisectionEC90
    if fa(c) - ystar < 0:
  File "/Users/*****/Documents/Python/NumericBisectionMethod/venv/twologisticfuncs.py", line 60, in f
    return c1 / (k1 + np.exp((-1*r1) * x))
KeyboardInterrupt

问题根源分析 #

1. ystar计算错误：
   原代码通过fa(x).max() - fa(x).min()计算函数值域，但x的采样区间[lo,up]可能不足以让logistic函数达到理论最值，导致ystar取值错误，二分法无法找到符合条件的解，陷入循环。
2. 二分法无终止边界：
   仅用abs(fa(c)-ystar) > 1e-9作为循环条件，当浮点精度无法满足该阈值时，会无限迭代；同时初始值c=1可能不在解的区间内，导致逻辑混乱。
3. 复合函数值域未正确计算：
   复合函数的最值与单个logistic函数不同，直接复用单个函数的二分逻辑会导致ystar偏差。

修复方案 #

1. 重构二分法函数（通用单logistic函数） #

直接使用logistic函数的理论最值计算ystar，添加迭代次数和区间长度双重终止条件：

def BisectionEC(c, k, r, a, b, target_frac):
    # 计算logistic函数的理论最值：max=c/k，min=0
    max_val = c / k
    min_val = 0
    ystar = target_frac * (max_val - min_val)
    
    def func(x):
        return c / (k + np.exp(-r * x))
    
    # 检查端点是否满足条件
    func_a = func(a)
    func_b = func(b)
    if abs(func_a - ystar) <= 1e-9:
        return a
    if abs(func_b - ystar) <= 1e-9:
        return b
    # 确保目标值在函数区间内
    if (func_a - ystar) * (func_b - ystar) >= 0:
        raise ValueError(f"目标值{ystar}不在函数区间[{func_a:.4f}, {func_b:.4f}]内，参数c={c},k={k},r={r}")
    
    iter_count = 0
    max_iter = 1000  # 最大迭代次数
    while (b - a) > 1e-10 and abs(func((a+b)/2) - ystar) > 1e-9 and iter_count < max_iter:
        mid = (a + b) / 2
        val = func(mid)
        if val < ystar:
            a = mid
        else:
            b = mid
        iter_count += 1
    
    if iter_count >= max_iter:
        print(f"警告：达到最大迭代次数，当前解{(a+b)/2:.6f}，函数值{func((a+b)/2):.6f}，目标{ystar:.6f}")
    return (a + b) / 2

2. 复合函数二分法单独实现 #

针对复合函数的最值特性，单独编写二分逻辑：

def BisectionEC_comp(c1, k1, r1, c2, k2, r2, a, b, target_frac):
    # 计算复合函数的理论最值
    g_max = c2 / k2
    comp_max = c1 / (k1 + np.exp(-r1 * g_max))
    g_min = 0
    comp_min = c1 / (k1 + np.exp(-r1 * g_min))
    
    ystar = target_frac * (comp_max - comp_min)
    
    def comp_func(x):
        return c1 / (k1 + np.exp(-r1 * (c2 / (k2 + np.exp(-r2 * x)))))
    
    # 检查端点
    comp_a = comp_func(a)
    comp_b = comp_func(b)
    if abs(comp_a - ystar) <= 1e-9:
        return a
    if abs(comp_b - ystar) <= 1e-9:
        return b
    if (comp_a - ystar) * (comp_b - ystar) >= 0:
        raise ValueError(f"目标值{ystar}不在复合函数区间[{comp_a:.4f}, {comp_b:.4f}]内")
    
    iter_count = 0
    max_iter = 1000
    while (b - a) > 1e-10 and abs(comp_func((a+b)/2) - ystar) > 1e-9 and iter_count < max_iter:
        mid = (a + b) / 2
        val = comp_func(mid)
        if val < ystar:
            a = mid
        else:
            b = mid
        iter_count += 1
    
    if iter_count >= max_iter:
        print(f"警告：复合函数达到最大迭代次数，当前解{(a+b)/2:.6f}，函数值{comp_func((a+b)/2):.6f}，目标{ystar:.6f}")
    return (a + b) / 2

3. 修改主逻辑中的调用代码 #

替换原有的EC计算部分：

# EC90 and EC10 for f(x), g(x) and f(g(x))
up = 20
lo = 0

EC90_1 = BisectionEC(c1, k1, r1, lo, up, 0.9)
EC10_1 = BisectionEC(c1, k1, r1, lo, up, 0.1)

EC90_2 = BisectionEC(c2, k2, r2, lo, up, 0.9)
EC10_2 = BisectionEC(c2, k2, r2, lo, up, 0.1)

EC90_3 = BisectionEC_comp(c1, k1, r1, c2, k2, r2, lo, up, 0.9)
EC10_3 = BisectionEC_comp(c1, k1, r1, c2, k2, r2, lo, up, 0.1)

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