Scipy.minimize：如何使用 epsilon 计算上确界

Question

我正在尝试确定以下函数的数值上界：

p(x) 图示例

p(x) 是一个概率函数，其中 p(0+) > 0，我进一步假设它是单调递减的。

问题是，在使用 epsilon 进行实验时，我没有得到一个好的值，我可以使用包含某些值的表来验证该值。 或者更确切地说，没有固定的 epsilon 可以为 p(x) 的所有实例提供最佳结果。

我的方法是使用

scipy.optimize.minimize

找到我的函数的最大值。

def density_upper_scaled(alpha, sigma):
    epsilon = 0.1
    constraint = {'type': 'ineq', 'fun': constraint_function, 'args': (alpha, sigma)}
    guess = np.array([1])
    result = sp.optimize.minimize(density_upper_function, guess, bounds=[(epsilon, None)], constraints=constraint)
    max_value = -result.fun  

    return (1.437 * 200**2 / max_value)
#, {'type': 'ineq', 'fun': lambda x: epsilon * x**2}

def density_upper_function(r):
    return - r**2 * 0.1


def constraint_function(r, alpha, sigma):
    return probability_function(r, alpha, sigma) - 0.1


def probability_function(distance, alpha, sigma):  
    ro0 = 200    # theoretical communication distance for sigma = 0
    if sigma != 0:
        p = 0.5 - 0.5 * math.erf((10 / math.sqrt(2)) * (alpha / sigma) * math.log10(distance / ro0))
    else:
        p = int(distance <= ro0)
    return p

我的问题是，如果 epsilon 太小，则最小化 -> inf。

然后我想自己用一个简单的脚本尝试一下，它有时表现更好，但对于 sigma/ alpha -> 0 来说表现更差，其想法如下。搜索最高的 x，以便我的 p(x) 仍然大于等于我的 epsilon，并且 epsilon r**2 也相同。

condition = probability_function(x, alpha, sigma)
    prior = -1
    epsilon = 0.2
    while condition >= epsilon and condition * x ** 2 >= epsilon:
        prior = condition
        condition = round(probability_function(x, alpha, sigma),5)

        x += 1
    print (x, condition)
    print(x,condition, 200**2 * 1.437 / (condition * x **2) )

我还尝试使用二阶导数来找到符号翻转的 x 或取其最大值的值，但它也与我需要的值不完全匹配。

有谁知道我的思维错误在哪里，或者有人知道另一种方法来实现我的目标吗？

Answer 1

让我们创建一个 MCVE 来解决您的问题：

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, special

我们从您的 Wolfram 链接中绘制概率函数（随意更改）：

def p(r, alpha=10., sigma=np.sqrt(2.), r0=200):  
    return 0.5 - 0.5 * special.erf(0.5 * (alpha / sigma) * np.log10(r / r0))

因此，您的模型是（其中减号是将上确界更改为下确界，因为我们将最小化而不是最大化）：

def model(x):
    # x = (epsilon, r)
    return - x[0] * np.power(x[1], 2.)

如果没有约束，解确实是无约束的。现在我们添加以下非线性约束：

def constraint(x):
    # x = (epsilon, r)
    return p(x[1]) - x[0]

nlc = optimize.NonlinearConstraint(constraint, 0., np.inf)

是时候绑定一切了：

solution = optimize.minimize(
     model, x0=[0.1, 1.],
     bounds=[(0., 1.), (0., np.inf)],
     constraints=[nlc]
)

返回一个绑定解：

     fun: -18497.79994916601
     jac: array([-64051.72680664,   -146.17871094])
 message: 'Optimization terminated successfully'
    nfev: 34
     nit: 14
    njev: 10
  status: 0
 success: True
       x: array([  0.28879471, 253.08442649])