Блог

Главная - Вопросы по программированию - Scipy минимизирует с ограничением для всего набора данных и всех значений, кроме нуля

Scipy минимизирует с ограничением для всего набора данных и всех значений, кроме нуля

#python #optimization #scipy #statistics #minimize

#python #оптимизация #сципи #Статистика #свести к минимуму

Вопрос:

Я пытаюсь получить параметр статистического уравнения для своих данных с помощью наиболее логарифмического метода. Уравнение, которое я хочу согласовать с ограничениями, выглядит следующим образом: уравнение и ограничения

Наилучшее соответствие достигается при минимуме моей функции правдоподобия l

Минимальный пример будет выглядеть так:

 import numpy as np

from numpy import log
from scipy.optimize import minimize

y =np.array([0.21863326, 0.19805154, 0.22953017, 0.21906749, 0.22067327,
       0.20638931, 0.17845443, 0.20008429, 0.21702199, 0.16334912,
       0.18480577, 0.17172182, 0.16495525, 0.15907978, 0.17029919,
       0.14020628, 0.16528562, 0.17141436, 0.14978351, 0.1329871 ,
       0.14036109, 0.15894933, 0.16783223, 0.17372222, 0.15986161,
       0.1654368 , 0.16348146, 0.15595923, 0.15192792, 0.12272897,
       0.17252942, 0.17164107, 0.16064716, 0.14564287, 0.14578649,
       0.14152733, 0.1354919 , 0.11175379, 0.1380746 , 0.12547517,
       0.15136653, 0.13984282, 0.18308302, 0.12271885, 0.15289988,
       0.13492309, 0.13499516, 0.13373476, 0.1034279 , 0.14278288,
       0.14574681, 0.11614764, 0.11256923, 0.14796558, 0.11459825,
       0.12417535, 0.15693744, 0.14159134, 0.11885544, 0.13164357,
       0.13445257, 0.13527885, 0.13472062, 0.12027512, 0.12072214,
       0.15361264, 0.12973932, 0.11003032, 0.13575847, 0.11980422,
       0.1187932 , 0.11152574, 0.14656588, 0.13885414, 0.13960315,
       0.12921241, 0.09522926, 0.14543513, 0.14980696, 0.11318417,
       0.10785905, 0.13858491, 0.11922434, 0.11760534, 0.12059705,
       0.12150726, 0.1184712 , 0.11084933, 0.10894509, 0.10107464,
       0.10258616, 0.1094653 , 0.095096  , 0.10059849, 0.10931144,
       0.11704954, 0.12639652, 0.13283708, 0.10203757, 0.10787873])

def l(para,args):
    xi,sigma = para
    y = args
    k = len(y)
    SUM = log(1 (xi*y)/sigma)
    return (-k*log(sigma)-(1 xi)/xi*np.sum(SUM))

def constrain1(x):
    xi,sigma = x
    return sigma
def constrain2(x):
    xi,sigma = x
    return -xi
def constrain3(x,*args):
    xi,sigma = x
    yi = args
    return 1 xi*yi/sigma
    

con1 = {'type': 'ineq', 'fun': constrain1}
con2 = {'type': 'ineq', 'fun': constrain2}
con3 = {'type': 'ineq', 'fun': constrain3,'args': y}
cons = [con1,con2,con3]

x0 =(-0.1,2)   
bound =((-np.inf,None),(None,np.inf))
res = minimize(l,x0,(y),bounds = bound,constraints=cons)

Первая проблема, с которой я столкнулся, — это реализовать ограничение 3, чтобы член внутри логарифма был больше нуля для каждой точки данных (1 (xi y_i) / sigma> 0). Я попытался добавить ограничения в цикл for

 con1 = {'type': 'ineq', 'fun': constrain1}
con2 = {'type': 'ineq', 'fun': constrain2}
cons = [con1,con2]
for yi in y:
    con3 = {'type': 'ineq', 'fun': constrain3,'args': yi}
    cons.append(con3)

но это тоже не сработало.

Моя вторая проблема заключается в том, что ограничения просто поддерживают >= или = . Есть ли способ реализовать просто> и особенно not equal (!=)? Я хочу xi ! = 0 и sigma> 0. Я очень удивлен, что ничего не нашел по этой теме, поскольку это должно быть стандартной проблемой в математике, или я недостаточно внимательно смотрел

Я нашел обходной путь для этого в mathlab, где я выполнил одну оптимизацию с xi<=0, а другую с xi>= 0 и выбрал лучшую. на сигму, похоже, не так сильно влияет ограничение, что она может быть равна нулю. Но я бы предпочел решение, в котором мне не нужно было бы переходить от одного решения к другому.

Комментарии:

1. Некоторые замечания: строгие неравенства не поддерживаются и теоретически не имеют особого смысла в большинстве оптимизационных фреймворков. Требуется явная формулировка на основе epsilon ( >= 0 eps ), и вы должны решить, является ли epsilon 10^-100 или 10^-2 . scipy.minimize в основном касается числовой оптимизации (и в лучшем случае: выпуклой), но x != 0 делит пространство на невыпуклый набор / дизъюнкцию, и вы не можете сделать это без методов ветвления / дискретной оптимизации. Упрощенно говоря: вы просите вещи, которые позволили бы нам легко решить все эти NP-сложные проблемы. Должен быть подвох.

2. Выбранный путь часто зависит от конкретного варианта использования, и некоторые знания, безусловно, помогают. Без контекста и чисто рассматривая ваши уравнения, я бы просто сказал, что x != 0 ограничение is — это хак, позволяющий избежать деления на ноль, и его никогда не следует обрабатывать как неравенство с разветвлением. В scipy optimize вы, вероятно, просто поймали бы этот случай и сделали что-то похожее на градиент. Если это сработает, это уже другая история.

3. Спасибо за объяснение @sascha . Вы правы. Казалось простым, что подвоха не будет. Я пытаюсь подогнать обобщенное распределение Парето к своим данным в контексте статистической оценки экстремальных значений, и ограничение `xi! = 0` действительно существует, поэтому вы не можете разделить на ноль. Это моя первая попытка использовать `scipy.minimize`, и мой математический опыт ограничен инженерной точкой зрения. Как я должен реализовать это в своем коде?

4. Я бы проигнорировал это и не использовал начальный запуск, где это значение не равно нулю. Шансы, что он получит ноль, довольно малы (и даже это может быть обработано в вашей функции с некоторым if; потеря плавности, надеюсь, не убьет решатель). Вам нужно немного поэкспериментировать. Я мог бы предположить, что вы столкнетесь с более чем одной числовой проблемой, включая некоторую необходимость настройки episilon для sigma. Кроме того, если он невыпуклый (не уверен), вы, вероятно, тоже захотите использовать несколько запусков.

Ответ №1:

 import numpy as np

from numpy import log
from scipy.optimize import minimize

y = np.array([0.21863326, 0.19805154, 0.22953017, 0.21906749, 0.22067327,
              0.20638931, 0.17845443, 0.20008429, 0.21702199, 0.16334912,
              0.18480577, 0.17172182, 0.16495525, 0.15907978, 0.17029919,
              0.14020628, 0.16528562, 0.17141436, 0.14978351, 0.1329871,
              0.14036109, 0.15894933, 0.16783223, 0.17372222, 0.15986161,
              0.1654368, 0.16348146, 0.15595923, 0.15192792, 0.12272897,
              0.17252942, 0.17164107, 0.16064716, 0.14564287, 0.14578649,
              0.14152733, 0.1354919, 0.11175379, 0.1380746, 0.12547517,
              0.15136653, 0.13984282, 0.18308302, 0.12271885, 0.15289988,
              0.13492309, 0.13499516, 0.13373476, 0.1034279, 0.14278288,
              0.14574681, 0.11614764, 0.11256923, 0.14796558, 0.11459825,
              0.12417535, 0.15693744, 0.14159134, 0.11885544, 0.13164357,
              0.13445257, 0.13527885, 0.13472062, 0.12027512, 0.12072214,
              0.15361264, 0.12973932, 0.11003032, 0.13575847, 0.11980422,
              0.1187932, 0.11152574, 0.14656588, 0.13885414, 0.13960315,
              0.12921241, 0.09522926, 0.14543513, 0.14980696, 0.11318417,
              0.10785905, 0.13858491, 0.11922434, 0.11760534, 0.12059705,
              0.12150726, 0.1184712, 0.11084933, 0.10894509, 0.10107464,
              0.10258616, 0.1094653, 0.095096, 0.10059849, 0.10931144,
              0.11704954, 0.12639652, 0.13283708, 0.10203757, 0.10787873])


def l(para, args):
    xi, sigma = para
    y = args
    k = len(y)
    # Define a list of functions to use np.sum() on
    SUM = [log(1   (xi * y_i) / sigma) for y_i in y]
    return -k * log(sigma) - (1   xi) / xi * np.sum(SUM)


def constrain1(x):
    (xi, sigma) = x
    return sigma


def constrain2(x):
    xi, sigma = x
    return -xi


# Use a class to prevent the python garbage collector to kick in
class Constraint3:
    def __init__(self, y_i):
        self.y_i = y_i

    def make_constraint(self, x):
        xi, sigma = x
        return 1   xi * self.y_i / sigma


con1 = {'type': 'ineq', 'fun': constrain1}
con2 = {'type': 'ineq', 'fun': constrain2}
constraints_list = [con1, con2]

for y_i in y:
    # Make each constraint function unique by instantiating the above class with the
    # new y_i value. The usage of individually defined function is not advised, because
    # some of them could get garbage collected, when there are too many items in the
    # array.
    constraints_list.append({'type': 'ineq', 'fun': Constraint3(y_i).make_constraint})

# Use proper ndarray-type
x0 = np.array([-0.1, 2])
bound = ((-np.inf, None), (None, np.inf))

res = minimize(l, x0, y, bounds=bound, constraints=constraints_list,
               options={'disp': True})
print("            Parameter values:")
print("              xi:   ", res.x[0])
print("              sigma:", res.x[1])

# Optimization terminated successfully    (Exit mode 0)
#         Current function value: -1739.2412627695364
#         Iterations: 24
#         Function evaluations: 82
#         Gradient evaluations: 24
#         Parameter values:
#           xi:    -1928.571453953235
#           sigma: 35762853.400437616
Метки: