In [1]:
SEED = sum(map(ord, "FOURIER_MSRA_2D"))
In [2]:
import sys
sys.path.append('../../..')
In [3]:
from lib.msra_analytical import AnalyticalLossFunctionAbs
In [4]:
import warnings
warnings.filterwarnings('error')
In [5]:
import numpy
import scipy
import scipy.stats
from scipy.integrate import quad, dblquad
In [6]:
def fourier_integral1d(cplx_integrand, j, x, eta):
def real_integrand(u):
return scipy.real(cplx_integrand(u, j, x, eta))
real_quad = quad(real_integrand, -numpy.inf, numpy.inf, epsrel=1e-4)[0]
return real_quad
In [7]:
def fourier_integral2d(cplx_integrand, j, k, x, y, eta):
def real_integrand(u, v):
return scipy.real(cplx_integrand(u, v, j, k, x, y, eta))
real_quad = dblquad(real_integrand,
-numpy.inf, numpy.inf,
lambda xx: -numpy.inf, lambda xx: numpy.inf,
epsrel=1e-4)[0]
return real_quad
In [8]:
class FourierGaussLossFunction(AnalyticalLossFunctionAbs):
def __init__(self, v_mu, v_sigma, alpha, c=None):
self.__v_mu = numpy.array(v_mu)
self.__v_sigma = numpy.array(v_sigma).reshape((len(v_mu), len(v_mu)))
self.__minus_1_half_over_pi = -0.5 / numpy.pi
self.__1_over_pi = 1. / numpy.pi
self.__1_over_2pi2 = 1 / (2 * numpy.pi)**2
self.__alpha = alpha
super(FourierGaussLossFunction, self).__init__(len(v_mu), c)
def moment_generating(self, t, i):
mu = self.__v_mu[i]
sigma2 = self.__v_sigma[i, i]
log_part = mu * t + 0.5 * sigma2 * t**2
return numpy.exp(log_part)
def moment_generating2D(self, vec_t, i, j):
if i > j:
i, j = j, i
vec_t = vec_t[::-1]
mu = self.__v_mu[[i, j]]
sigma2 = [[self.__v_sigma[i, i], self.__v_sigma[i, j]],
[self.__v_sigma[j, i], self.__v_sigma[j, j]]]
log_part = numpy.dot(mu, vec_t)
log_part += 0.5 * numpy.dot(vec_t, numpy.dot(sigma2, vec_t))
return numpy.exp(log_part)
def e(self, m):
return (self.__v_mu - m).sum()
def g_fourier_integrand(self, u, j, m_j, eta):
i = 1j
eta_m_iu = eta - i * u
res = numpy.exp(-eta_m_iu * m_j)
res *= self.moment_generating(eta_m_iu, j)
res /= i * (u + i * eta)**3
return res
def g(self, i, m_i):
continue_bool = True
eta = 1.5 * numpy.random.rand()
while continue_bool:
try:
integral = fourier_integral1d(self.g_fourier_integrand, i, m_i, eta)
continue_bool = False
return self.__1_over_pi * integral
except scipy.integrate.IntegrationWarning, e:
#print "g not converging for x = %s, alpha = %s" % (m_i, eta)
eta = 1.5 * numpy.random.rand()
def f_fourier_integrand(self, u, j, m_j, eta):
i = 1j
eta_m_iu = eta - i * u
res = numpy.exp(-eta_m_iu * m_j)
res *= self.moment_generating(eta_m_iu, j)
res /= (u + i * eta)**2
return res
def f(self, i, m_i):
continue_bool = True
eta = 1.5 * numpy.random.rand()
while continue_bool:
try:
integral = fourier_integral1d(self.f_fourier_integrand, i, m_i, eta)
continue_bool = False
return self.__minus_1_half_over_pi * integral
except scipy.integrate.IntegrationWarning, e:
#print "f not converging for x = %s, alpha = %s" % (m_i, eta)
eta = 1.5 * numpy.random.rand()
def h_fourier_integrand(self, u, v, j, k, x, y, eta):
i = 1j
eta_m_iu = eta - i * numpy.array([u, v])
res = numpy.exp(numpy.dot(-eta_m_iu, [x, y]))
res *= self.moment_generating2D(eta_m_iu, j, k)
res /= (u + i*eta[0])**2 * (v + i*eta[1])**2
return res
def h(self, i, j, m_i, m_j):
continue_bool = True
eta = 1.5 * numpy.random.rand(2)
while continue_bool:
try:
integral = fourier_integral2d(self.h_fourier_integrand, i, j, m_i, m_j, eta)
continue_bool = False
return self.__1_over_2pi2 * integral
except scipy.integrate.IntegrationWarning, e:
eta = 1.5 * numpy.random.rand(2)
def l_fourier_integrand(self, u, v, j, k, x, y, eta):
i = 1j
eta_m_iu = eta - i * numpy.array([u, v])
res = numpy.exp(numpy.dot(-eta_m_iu, [x, y]))
res *= self.moment_generating2D(eta_m_iu, j, k)
res /= (u + i*eta[0])**2 * (-eta_m_iu[1])
return res
def l(self, i, j, m_i, m_j):
continue_bool = True
eta = 1.5 * numpy.random.rand(2)
while continue_bool:
try:
integral = fourier_integral2d(self.l_fourier_integrand, i, j, m_i, m_j, eta)
continue_bool = False
return self.__1_over_2pi2 * integral
except scipy.integrate.IntegrationWarning, e:
eta = 1.5 * numpy.random.rand(2)
def shortfall_risk(self, m=None):
m = self._check_argument(m)
sum_e = self.e(m)
sum_g, sum_h = 0., 0.
for i, m_i in enumerate(m):
sum_g += self.g(i, m_i)
for j, m_j in enumerate(m):
if j > i:
sum_h += self.h(i, j, m_i, m_j)
return sum_e + 0.5 * sum_g + self.__alpha * sum_h
def shortfall_risk_jac(self, m):
m = self._check_argument(m)
res = []
for i, m_i in enumerate(m):
partial_der = 1 + self.f(i, m_i)
for j, m_j in enumerate(m):
if i != j:
partial_der += self.__alpha * self.l(j, i, m_j, m_i)
res.append(partial_der)
return numpy.array(res)
In [9]:
rho = 0.2
mu = [0., 0., 0.]
sigma = [[0.5, 0.5 * rho, 0.], [0.5 * rho, 0.5, 0.], [0., 0., 0.6]]
c = 1.
loss = FourierGaussLossFunction(mu, sigma, 1., c)
In [10]:
from scipy.optimize import minimize
maxiter = 3500
In [11]:
guess = 0. * numpy.ones(loss.dim)
In [12]:
cons = ({'type': 'ineq',
'fun' : lambda x: loss.ineq_constraint(x),
'jac' : lambda x: loss.ineq_constraint_jac(x)})
In [13]:
numpy.random.seed(SEED)
res = minimize(loss.objective, guess,
jac=loss.objective_jac,
constraints=cons,
method='SLSQP',
options={'maxiter': maxiter})
In [14]:
print res
print
print loss.ineq_constraint(res.x)
print numpy.mean(res.x[:-2])
In [ ]: