Optimisation de la variance des variables bornées

Question

Soit x_1, ..., x_i, ..., x_p p p nombres réels tels que 0 < = x_i < = b pour tout i. En d'autres termes, chaque x_i peut prendre n'importe quelle valeur entre 0 et b. Je voudrais trouver les valeurs de {x_i} qui maximisent la variance entre elles.

Avez-vous des conseils? J'aimerais utiliser ce résultat pour mon exemple de code. Ou, cette question n'est-elle pas bien définie? Au début, j'ai pensé à quelque chose comme x_1 = 0, x_2 = ... = x_p = b, mais j'ai ensuite trouvé que cela ne maximise pas la variance lorsque p est un peu grand.

Merci

Answer 1

Après les commentaires, je l'ai fait quelques essais sur un numérique prouvent votre problème. Il y a encore du travail à faire, mais j'espère que cela vous mettra sur la bonne voie. En outre, j'ai utilisé python, et je n'ai aucune idée si cela est bon pour vous ou non. Vous pouvez certainement trouver des façons équivalentes de le faire en matlab et R. J'utilise la propriété bien connue de la variance = E [X^2] - E [X]^2, pour faciliter les dérivées. (Si vous avez des doutes, vérifiez wiki).

Le package pythonscipy.optimize a une méthode minimize pour minimiser numériquement une fonction. Vous pouvez sélectionner un algorithme pour résoudre le problème. Je ne suis pas si familier avec les algorithmes possibles et je cherchais la descente en dégradé uni bien connue (enfin, j'espère que vous le savez), et je pense que celle qui est fermée pourrait être SLSQP, mais honnêtement je ne suis pas 100% sûr sur les détails. Enfin, je ne m'assurais pas que la fonction que vous minimisez est convexe, ou je me suis rendu compte si elle avait des minima locaux, mais les résultats sont corrects.

Je vous donne le code en python ci-dessous, dans le cas où il est utile, mais le fond est que je vous propose:

• Choisissez une langue/package que vous êtes familier avec
• Choisissez une algorithme pour l'optimisation
• ce serait bien de prouver que la fonction est convexe (de sorte que la solution converge)
• Définissez les paramètres que vous voulez faire le prouver

Code ci-dessous. J'espère que cela aide.

Je ne vais pas publier l'algèbre pour les dérivés, j'espère que vous pouvez les fabriquer vous-même. Et vous devez prendre en compte le fait que vous maximisez et ne minimisez pas, donc vous devez multiplier par -1, comme expliqué très clairement j'espère que here (cherchez "maximiser").

Setup,

In [1]: 

from scipy.optimize import minimize 
import numpy as np 

La fonction vous maximisez, qui est la variance (rappelez-vous l'affaire E [X^2] - E [X]^2, et -1),

In [86]: 

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

le dérivé de cette fonction pour chacun des xi du vecteur x, (je l'espère, vous pouvez obtenir et à dériver le même résultat),

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) 

En fait, j'ai fait pas mal d'erreurs en écrivant cette fonction, à la fois dans l'implémentation dérivée et python. Mais il y a une astuce qui aide beaucoup, qui est de vérifier la dérivée d'une manière numérique, en ajoutant et en soustrayant un petit epsilon dans chaque dimension et en calculant la pente de la courbe see wiki. Cela la fonction qui se rapproche du dérivé,

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 

Et puis je l'ai vérifié func_deriv_approx contre func_deriv pour un tas de valeurs.

Et la minimisation elle-même. Si j'initialiser les valeurs à la solution que nous soupçonnons est juste, il fonctionne bien, il ne itère une fois et donne les résultats escomptés,

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.] 

(Notez que vous pouvez utiliser la longueur que vous vouliez, depuis func et func_deriv sont écrit d'une manière qui accepte n'importe quelle longueur).

Vous pouvez initialiser au hasard comme ça,

In [81]: 

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

In [82]: 

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

Et puis la maximisation est,

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] 

Ou enfin pour la longueur = 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]