In [1]:

    

import sys

sys.path.append('../../..')


    

In [2]:

    

import numpy
import scipy.stats

from lib.msra_loss import MSRALossFunctionAbs


    

In [3]:

    

class MCLossFunction(MSRALossFunctionAbs):
    def __init__(self, distrib, alpha, c=None):
        self.__alpha = alpha
        super(MCLossFunction, self).__init__(distrib, c)
        
    def shortfall_risk(self, m=None):
        m = self._check_argument(m)
        x_minus_m = numpy.subtract(self.x, m)
        
        mean_sum_ = numpy.mean(x_minus_m.sum(axis=1))
        
        pos_part = numpy.maximum(x_minus_m, 0.)
        pos_part_squared = numpy.square(pos_part)
        mean_sum_2_ = numpy.mean(pos_part_squared.sum(axis=1))
        
        to_add = 0.
        for i in xrange(self.dim):
            for j in xrange(i + 1, self.dim):
                to_add += numpy.mean(numpy.multiply(pos_part[:, i], pos_part[:, j]))
                        
        return mean_sum_ + 0.5 * mean_sum_2_ + self.__alpha * to_add

    def shortfall_risk_jac(self, m):
        m = self._check_argument(m)
        x_minus_m = numpy.subtract(self.x, m)
        
        pos_part = numpy.maximum(x_minus_m, 0.)
        mean_pos_part = numpy.mean(pos_part, axis=0)
        
        dbl = []
        for i in xrange(self.dim):            
            indic_i = numpy.sign(pos_part[:, i])
            tmp = 0.
            for j in xrange(self.dim):
                if i != j:
                    tmp += numpy.mean(numpy.multiply(indic_i, pos_part[:, j]))
                
            dbl.append(self.__alpha * tmp)
        
        return mean_pos_part + 1. + dbl


    

In [4]:

    

M = 100000

rho = -0.5

mu = [0., 0., 0.]
sigma = [[0.5, 0.5 * rho, 0.], [0.5 * rho, 0.5, 0.], [0., 0., 0.6]]

rv = scipy.stats.multivariate_normal(mean=mu, cov=sigma, allow_singular=True)
X = rv.rvs(size=M)

c = 1.

loss = MCLossFunction(X, 1., 1.)


    

In [5]:

    

from scipy.optimize import minimize

maxiter = 3500


    

In [6]:

    

guess = 0. * numpy.ones(loss.dim)


    

In [7]:

    

cons = ({'type': 'ineq',
         'fun' : lambda x: loss.ineq_constraint(x),
         'jac' : lambda x: loss.ineq_constraint_jac(x)})


    

In [8]:

    

res = minimize(loss.objective, guess, 
               jac=loss.objective_jac, 
               constraints=cons, 
               method='SLSQP',
               options={'maxiter': maxiter, 'disp': True})


    

    




Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.255311639591
            Iterations: 5
            Function evaluations: 5
            Gradient evaluations: 5

In [9]:

    

print res
print
print loss.ineq_constraint(res.x)
print numpy.mean(res.x[0:1])


    

    




     fun: -0.25531163959115105
     jac: array([ 1.,  1.,  1.,  0.])
 message: 'Optimization terminated successfully.'
    nfev: 5
     nit: 5
    njev: 5
  status: 0
 success: True
       x: array([-0.13662206, -0.13510257,  0.01641299])

-7.69532349842e-07
-0.136622057326