2015-09-11 31 views
3

我有一個數據集,我想通過最小平方誤差方法找到混合高斯模型。Scipy&Optimize:最小化示例,如何添加約束?

的代碼是這樣的:

from sklearn.neighbors import KernelDensity 
kde = KernelDensity().fit(sample) 

def gaussian_2d(x,y,meanx,meany,sigx,sigy,rho): 
    # rho <= 1 
    part1 = 1/(2*np.pi*sigx*sigy*sqrt(1-0.5**2)) 
    part2 = -1/2*(1-rho**2) 
    part3 = (((x-meanx)/sigx)**2-2*rho*(x-meanx)*(y-meany)/(sigx*sigy)+((y-meany)/sigy)**2) 
    return part1*exp(part2*part3) 

def square_error(f1,f2, u1,v1,sigu1,sigv1,rho1, u2,v2,sigu2,sigv2,rho2, u3,v3,sigu3,sigv3,rho3): 
    # 1. Generate Mixed Gaussian Model 
    def gaussian1(x,y): 
     return gaussian_2d(x,y,u1,v1,sigu1,sigv1,rho1) 
    def gaussian2(x,y): 
     return gaussian_2d(x,y,u2,v2,sigu2,sigv2,rho2) 
    def gaussian3(x,y): 
     return gaussian_2d(x,y,u3,v3,sigu3,sigv3,rho3) 
    mixed_model = f1*gaussian1(x,y)+f2*gaussian2(x,y)+(1-f1-f2)*gaussian3(x,y) 
    # 2. Calculate the sum of square error 
    sum_error = 0 
    for point in sample: 
     error = (exp(mixed_model(point)) - exp(kde.score(point)))**2 
     sum_error += error 
    return sum_error 

# How can I add constraints: 
# f1+f2 <= 1 
# rho1,2,3 <= 1 
result = sp.optimize.minimize(square_error) 

但我不知道如何在minimize方法添加約束。 http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize中的示例很難理解。

UPDATE: 這是我結束了,

result = sp.optimize.minimize(
    square_error, 
    x0 = [0.2,0.5, 
      1,1,1,1,0.3, 
      1,1,1,1,0.3, 
      1,1,1,1,0.3,], 
    bounds = [(0., 1.),(0., 1.), 
       (None, None),(None, None),(0., None),(0., None),(0., 1.), 
       (None, None),(None, None),(0., None),(0., None),(0., 1.), 
       (None, None),(None, None),(0., None),(0., None),(0., 1.),]) 

但它給了我TypeError: square_error() takes exactly 17 arguments (1 given),有什麼問題呢?

+0

可能重複的[scipy最小化約束](http://stackoverflow.com/questions/20075714/scipy-minimize-with-constraints) –

回答

1

如果解算器支持它們,則只能添加添加邊界,因此僅適用於method='L-BFGS-B'TNCSLSQP

通過一系列(min, max)元組來傳遞邊界,這些元組的長度對應於您的參數數量。擬合3個參數的示例如下:

result = sp.optimize.minimize(
         square_error, 
         method='L-BFGS-B', 
         bounds=[(0., 5.), (None, 1.e4), (None, None)]) 

這裏,None對應於無界限。 恐怕您的示例中f1+f2 <= 1等參數組合的限制在的bounds框架內是不可能的。

但是,如果違反了您的界限,您可以簡單地在您的成本函數中返回np.inf。我不知道這樣做的穩定性,雖然:

def square_error(f1,f2, other_args): 
    if f1+f2 <= 1: 
     return np.inf 
    # rest of the cost function 

此外,我建議使用,而不是從頭開始創建他們的multivariate Gaussians Python實現。這將加快你的裝修,幫助避免錯誤並提高可讀性。

+0

'我恐怕約束組合的參數,如f1 +在你的例子中f2 <= 1是不可能的。「這是不正確的。有一些解算器特別適用於支持它們的約束優化。 – cel

+0

你能提供一個鏈接嗎?我指的是包含在'scipy.minimize'中的求解器。 –

+1

只需看一下'constraints'關鍵字即可。 OP正試圖定義約束而不是簡單的界限。 – cel