【问题标题】:Maximizing variance of bounded variables最大化有界变量的方差
【发布时间】:2015-09-05 20:44:48
【问题描述】:

令 x_1, ..., x_i,..., x_p 是 p 个实数,使得 0

你有什么提示吗? 我想将此结果用于我的示例代码。 或者,这个问题不是很好定义的吗?

一开始我想到了 x_1=0, x_2=...=x_p=b 之类的东西,但后来我发现当 p 有点大时,这并不能最大化方差。

谢谢

【问题讨论】:

  • 不应该是一半的值等于0,另一半等于b吗?
  • 谢谢,这对我来说听起来很合理。你知道如何证明这一点吗?
  • 起来!我的代数有点生疏......(此外,我必须学习如何使用符号来写答案!)刚刚拿了一张纸,并且只能证明,对于偶数 p,方差对于像上面这样选择的值将是 (1/2)*b^2。也许用数值方法......但这真的很奇怪!你真的需要证明吗?
  • 好吧,我不想给您带来太多麻烦,但是如果您能概述一下您的流程,或者将我指向一个解释这一点的链接,那将非常有帮助。我真的很感谢你的帮助!谢谢。

标签: statistics variance maximization


【解决方案1】:

在 cmets 之后,我对您的问题进行了一些数值证明的试验。还有一些工作要做,但我希望它能让你走上正轨。此外,我使用过python,我不知道这是否适合您。您肯定可以在 matlabR 中找到等效的方法。

我使用方差 = E[X^2] - E[X]^2 的众所周知的性质来简化导数。 (如有疑问,请查看wiki)。

pythonscipy.optimize 有一个方法minimize 可以在数值上最小化一个函数。您可以选择解决问题的算法;我对可能的算法不太熟悉,我一直在寻找众所周知的普通梯度下降(好吧,至少我希望你知道),我认为封闭的可能是SLSQP,但老实说我'我不能 100% 确定细节。

最后,我不确定您要最小化的函数是否是凸函数,或者确定它是否具有局部最小值,但结果看起来不错。

我给你下面的python代码,以防它有用,但最重要的是我建议你:

  • 选择您熟悉的语言/包
  • 选择优化算法
  • 最好证明该函数是凸函数(以便解收敛)
  • 设置要验证的参数

代码如下。希望对您有所帮助。

我不打算发布导数的代数,希望您可以自己制作。而且你必须考虑到你是在最大化而不是最小化,所以你必须乘以 -1,正如我希望很清楚的解释 here(寻找“最大化”)。

设置,

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)

对于向量x 的每个xi 的该函数的导数,(我希望您可以推导并得到相同的结果),

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 实现中。但是有一个技巧很有帮助,那就是以数字方式检查导数,通过在每个维度上加减一个小 epsilon 并计算曲线的斜率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_approxfunc_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.]

(请注意,您可以使用所需的长度,因为 funcfunc_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]

【讨论】:

  • 如果您已经阅读了我的回答;我终于找到了问题所在,因此我将对其进行大量编辑,而不尊重我的原始文本。请过几分钟再检查。
  • 这真的很有帮助,非常感谢!!是的,我了解 Python,并且基本上复制了您的想法。也感谢您提供有关衍生品的信息。现在我想我会再次研究代数。
  • 现在我做到了!抱歉我不是很习惯这里的系统,再次非常感谢!
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2018-08-08
  • 2022-01-18
  • 2011-10-24
  • 1970-01-01
  • 2022-01-16
  • 2021-12-13
  • 2012-10-22
相关资源
最近更新 更多