Максимизация дисперсии ограниченных переменных

Question

Пусть x_1, ..., x_i, ..., x_p являются p действительными числами, такими, что 0 < = x_i < = b для всех i. То есть каждый x_i может принимать любое значение между 0 и b. Я хотел бы найти значения {x_i}, которые максимизируют дисперсию между ними.

Есть ли у вас какие-либо намеки? Я бы хотел использовать этот результат для моего примера кода. Или, не этот вопрос четко определен?

Сначала я подумал о чем-то вроде x_1 = 0, x_2 = ... = x_p = b, но тогда я обнаружил, что это не максимизирует дисперсию, когда p немного велико.

Благодаря

Answer 1

После замечаний, я сделал некоторые испытания на численном доказать для вашей проблемы. Есть еще какая-то работа, но я надеюсь, что это положит вас на правильный путь. Кроме того, я использовал python, и я понятия не имею, нормально ли это для вас или нет. Вы можете найти эквивалентные способы сделать это в matlab и R.

Я использую известное свойство дисперсии = E [X^2] - E [X]^2, чтобы упростить производные. (Если у вас есть сомнения, проверьте wiki).

python упаковка scipy.optimize имеет способ minimize, чтобы минимизировать числовую функцию. Вы можете выбрать алгоритм для решения проблемы; Я не очень хорошо знаком с возможными алгоритмами, и я искал хорошо известный градиентный спуск (ну, по крайней мере, я надеюсь, вы это знаете), и я думаю, что закрытым может быть SLSQP, но, честно говоря, я не 100% уверены в деталях.

И, наконец, я не удостоверился, что функция, которую вы минимизируете, является выпуклой или вычисляется, имеет ли она локальные минимумы, но результаты выглядят нормально.

Я даю вам код в питона ниже, в случае, если это полезно, но Bottomline в том, что я хотел бы предложить вам:

• Выберите язык/пакет вы знакомы с
• Выберите один алгоритм оптимизации
• было бы хорошо, чтобы доказать, что функция выпукла (так что решение сходится)
• Установите параметры, которые вы хотите сделать доказать

Код ниже. Надеюсь, поможет.

Я не собираюсь публиковать алгебру для производных, надеюсь, вы сможете сделать их сами. И вы должны принять во внимание, что вы максимизируете и не минимизируете, поэтому вам нужно умножить на -1, как объяснялось, я надеюсь, совершенно ясно here (ищите «максимизацию»).

Setup,

In [1]: 

from scipy.optimize import minimize 
import numpy as np 

Функция вы максимально, то есть дисперсия (помните трюк E [X^2] - E [X]^2, и -1),

In [86]: 

def func(x): 
    return (-1) * (np.mean([xi**2 for xi in x]) - np.mean(x)**2) 

производная этой функции для каждого из xi вектора x, (я надеюсь, что вы можете получить и порождают к тому же результату),

In [87]: 

def func_deriv(x): 
    n = len(x) 
    l = [] 
    for i in range(n): 
     res = (2 * x[i]/n) - ((2/(n**2)) * (x[i] + sum([x[j] for j in range(n) if j != i]))) 
     l += [(-1) * res] 
    return np.array(l) 

На самом деле, я допустил немало ошибок при написании этой функции как в производной, так и в реализации python. Но есть трюк, который помогает много, чтобы проверить производную по-цифровому, путем добавления и вычитания небольшого эпсилона в каждом измерении и вычисления наклона кривой see wiki. Это будет функция, которая приближается к производной,

In [72]: 

def func_deriv_approx(x, epsilon=0.00001): 
    l = [] 
    for i in range(len(x)): 
     x_plus = [x[j]+((j == i)*epsilon) for j in range(len(x))] 
     x_minus = [x[j]-((j == i)*epsilon) for j in range(len(x))] 
     res = (-1) * (func(x_plus) - func(x_minus))/(2*epsilon) 
     l += [res] 
    return l 

А потом я проверил func_deriv_approx против func_deriv для связки значений.

И минимизация себя. Если я инициализировать значение к решению мы подозреваем правильно, она работает хорошо, это только итерацию один раз и дает ожидаемый результат,

In [99]: 

res = minimize(func, [0, 0, 10, 10], jac=func_deriv, bounds=[(0,10) for i in range(4)], 
       method='SLSQP', options={'disp': True}) 
Optimization terminated successfully. (Exit mode 0) 
      Current function value: -25.0 
      Iterations: 1 
      Function evaluations: 1 
      Gradient evaluations: 1 

In [100]: 

print(res.x) 
[ 0. 0. 10. 10.] 

(Обратите внимание, что вы можете использовать длину вы хотели, так как func и func_deriv являются написанные таким образом, чтобы они принимали любую длину).

Вы можете инициализировать случайно, как это,

In [81]: 

import random 
xinit = [random.randint(0, 10) for i in range(4)] 

In [82]: 

xinit 
Out[82]: 
[1, 2, 8, 7] 

И тогда максимизация,

In [83]: 

res = minimize(func, xinit, jac=func_deriv, bounds=[(0,10) for i in range(4)], 
       method='SLSQP', options={'disp': True}) 
Optimization terminated successfully. (Exit mode 0) 
      Current function value: -25.0 
      Iterations: 3 
      Function evaluations: 3 
      Gradient evaluations: 3 
In [84]: 

print(res.x) 
[ 1.27087156e-13 1.13797860e-13 1.00000000e+01 1.00000000e+01] 

Или, наконец, для длины = 100,

In [85]: 

import random 
xinit = [random.randint(0, 10) for i in range(100)] 

In [91]: 

res = minimize(func, xinit, jac=func_deriv, bounds=[(0,10) for i in range(100)], 
       method='SLSQP', options={'disp': True}) 
Optimization terminated successfully. (Exit mode 0) 
      Current function value: -24.91 
      Iterations: 23 
      Function evaluations: 22 
      Gradient evaluations: 22 
In [92]: 

print(res.x) 
[ 2.49143492e-16 1.00000000e+01 1.00000000e+01 -2.22962789e-16 
    -3.67692105e-17 1.00000000e+01 -8.83129256e-17 1.00000000e+01 
    7.41356521e-17 3.45804774e-17 -8.88402036e-17 1.31576404e-16 
    1.00000000e+01 1.00000000e+01 1.00000000e+01 1.00000000e+01 
    -3.81854094e-17 1.00000000e+01 1.25586928e-16 1.09703896e-16 
    -5.13701064e-17 9.47426071e-17 1.00000000e+01 1.00000000e+01 
    2.06912944e-17 1.00000000e+01 1.00000000e+01 1.00000000e+01 
    -5.95921560e-17 1.00000000e+01 1.94905365e-16 1.00000000e+01 
    -1.17250430e-16 1.32482359e-16 4.42735651e-17 1.00000000e+01 
    -2.07352528e-18 6.31602823e-17 -1.20809001e-17 1.00000000e+01 
    8.82956806e-17 1.00000000e+01 1.00000000e+01 1.00000000e+01 
    1.00000000e+01 1.00000000e+01 3.29717355e-16 1.00000000e+01 
    1.00000000e+01 1.00000000e+01 1.00000000e+01 1.00000000e+01 
    1.43180544e-16 1.00000000e+01 1.00000000e+01 1.00000000e+01 
    1.00000000e+01 1.00000000e+01 2.31039883e-17 1.06524134e-16 
    1.00000000e+01 1.00000000e+01 1.00000000e+01 1.00000000e+01 
    1.77002357e-16 1.52683194e-16 7.31516095e-17 1.00000000e+01 
    1.00000000e+01 3.07596508e-17 1.17683979e-16 -6.31665821e-17 
    1.00000000e+01 2.04530928e-16 1.00276075e-16 -1.20572493e-17 
    -3.84144993e-17 6.74420338e-17 1.00000000e+01 1.00000000e+01 
    -9.66066818e-17 1.00000000e+01 7.47080743e-17 4.82924982e-17 
    1.00000000e+01 -9.42773478e-17 1.00000000e+01 1.00000000e+01 
    1.00000000e+01 1.00000000e+01 1.00000000e+01 5.01810185e-17 
    -1.75162038e-17 1.00000000e+01 6.00111991e-17 1.00000000e+01 
    1.00000000e+01 7.62548028e-17 -6.90706135e-17 1.00000000e+01]