Python二分法求EC10/EC90报only size-1 arrays类型错误

报错根因 #

触发TypeError: only size-1 arrays can be converted to Python scalars的直接原因是二分法函数内的语句x = float(np.linspace(a,b,1000))逻辑错误：
np.linspace(a,b,1000)返回的是包含1000个元素的一维numpy数组，而Python内置float()仅支持将单个尺寸为1的数值/数组转换为浮点标量，无法处理多元素数组。

除了这个直接报错，代码中还存在3处会导致计算结果完全错误的逻辑问题，一并修复即可正常运行：

• EC90目标值计算错误：当前BisectionEC90函数中ystar取的是10%幅值，和EC10的目标值完全一致，不符合EC90的定义
• 二分法循环条件写反：当前循环仅当函数值和目标值误差小于阈值时才执行，初始状态误差大于阈值时循环根本不会启动，直接返回初始值c=1
• 绘图变量赋值错误：提前将x赋值为单个数值1，无法绘制连续函数曲线

具体修复操作 #

• 移除两个二分法函数中np.linspace外层多余的float()转换，直接使用生成的数组计算函数在区间内的最大、最小值
• 修正BisectionEC90的目标值计算公式，将系数从0.1改为0.9
• 反转二分法循环的判断条件，当函数值和目标值的误差大于精度阈值（1e-9）时持续迭代缩小区间
• 把绘图用的x变量改回linspace生成的连续数组，不要赋值为单个标量1
• 优化初始迭代点赋值：将二分法初始中点设为搜索区间中点，替代固定值1，避免初始点不在搜索区间导致的收敛失败；补充EC值计算的最小值偏移，兼容函数最小值不为0的场景

修复后可运行代码 #

# Imports
import numpy as np
from numpy import log as ln
from matplotlib import pyplot as plt
from random import randint
from mpmath import *
import sympy as sym


# Params
u = 10
c1 = randint(1,u)
n1 = randint(1,u)
k1 = randint(1,u)

c2 = randint(1,u)
n2 = randint(1,u)
k2 = randint(1,u)

print(f"c1 = {c1}")
print(f"n1 = {n1}")
print(f"k1 = {k1}")

print('{}(x^{})/({}+x^{})'.format(c1,n1,k1,n1))

print(f"c2 = {c2}")
print(f"n2 = {n2}")
print(f"k2 = {k2}")

print('{}(x^{})/({}+x^{})'.format(c2,n2,k2,n2))

mp.dps = 15
mp.pretty = False

# 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*(x**n1)/(k1+x**n1)
def g(x):
    return c2*(x**n2)/(k2+x**n2)
def c(x):
    return ((c1*x**(n2*n1))*(c2**n1))/((k2*(k2+(x**n2))**n1)+((c2**n1)*(x**(n2*n1))))

# EC finder
def BisectionEC10(fa,a,b):
    c = (a + b)/2
    x = np.linspace(a,b,1000)
    y_min = fa(x).min()
    y_range = fa(x).max() - y_min
    ystar = 0.10 * y_range + y_min
    while abs(fa(c)-ystar) > 1e-9:
        c=(a+b)/2
        if fa(c)-ystar<0:
            a=c
        elif fa(c)-ystar>0:
            b=c
    return c


def BisectionEC90(fa,a,b):
    c = (a+b)/2
    x = np.linspace(a,b,1000)
    y_min = fa(x).min()
    y_range = fa(x).max() - y_min
    ystar= 0.9 * y_range + y_min
    while abs(fa(c)-ystar) > 1e-9:
        c=(a+b)/2
        if fa(c)-ystar<0:
            a=c
        elif fa(c)-ystar>0:
            b=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)

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(c, lo, up)
EC10_3 = BisectionEC10(c, lo, up)

# Hill Coefficient for f(x) and g(x)
H_1 = ln(81)/(ln(EC90_1/EC10_1))
H_2 = ln(81)/(ln(EC90_2/EC10_2))
H_3 = ln(81)/(ln(EC90_3/EC10_3))
print(f"Hill's Coefficient for f(x) = {H_1}")
print(f"Hill's Coefficient for g(x) = {H_2}")
print(f"Hill's Coefficient for f(g(x)) = {H_3}")

# Plotting
plt.plot(x, f(x), color = 'red', label='f(x)')
plt.plot(x, g(x), color = 'blue', label='g(x)')
plt.plot(x, c(x), color = 'green', label='f(g(x))')
plt.legend()
plt.show()

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