Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 49 additions & 50 deletions pykliep.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,35 +3,35 @@

class DensityRatioEstimator:
"""
Class to accomplish direct density estimation implementing the original KLIEP
Class to accomplish direct density estimation implementing the original KLIEP
algorithm from Direct Importance Estimation with Model Selection
and Its Application to Covariate Shift Adaptation by Sugiyama et al.
The training set is distributed via
and Its Application to Covariate Shift Adaptation by Sugiyama et al.

The training set is distributed via
train ~ p(x)
and the test set is distributed via
and the test set is distributed via
test ~ q(x).

The KLIEP algorithm and its variants approximate w(x) = q(x) / p(x) directly. The predict function returns the
estimate of w(x). The function w(x) can serve as sample weights for the training set during
training to modify the expectation function that the model's loss function is optimized via,
i.e.

E_{x ~ w(x)p(x)} loss(x) = E_{x ~ q(x)} loss(x).
Usage :
The fit method is used to run the KLIEP algorithm using LCV and returns value of J
trained on the entire training/test set with the best sigma found.

Usage :
The fit method is used to run the KLIEP algorithm using LCV and returns value of J
trained on the entire training/test set with the best sigma found.

Use the predict method on the training set to determine the sample weights from the KLIEP algorithm.
"""
def __init__(self, max_iter=5000, num_params=[.1,.2], epsilon=1e-4, cv=3, sigmas=[.01,.1,.25,.5,.75,1], random_state=None, verbose=0):
"""

def __init__(self, max_iter=5000, num_params=[.1,.2], epsilon=1e-4, cv=7, sigmas=[.01,.1,.25,.5,.75,1], random_state=None, verbose=0):
"""
Direct density estimation using an inner LCV loop to estimate the proper model. Can be used with sklearn
cross validation methods with or without storing the inner CV. To use a standard grid search.


max_iter : Number of iterations to perform
num_params : List of number of test set vectors used to construct the approximation for inner LCV.
Must be a float. Original paper used 10%, i.e. =.1
Expand All @@ -45,29 +45,29 @@ def __init__(self, max_iter=5000, num_params=[.1,.2], epsilon=1e-4, cv=3, sigmas
self.sigmas = sigmas
self.cv = cv
self.random_state = 0

def fit(self, X_train, X_test, alpha_0=None):
""" Uses cross validation to select sigma as in the original paper (LCV).
In a break from sklearn convention, y=X_test.
The parameter cv corresponds to R in the original paper.
Once found, the best sigma is used to train on the full set."""

# LCV loop, shuffle a copy in place for performance.
cv = self.cv
chunk = int(X_test.shape[0]/float(cv))
if self.random_state is not None:
np.random.seed(self.random_state)
X_test_shuffled = X_test.copy()
np.random.shuffle(X_test_shuffled)
np.random.permutation(X_test_shuffled)

j_scores = {}

if type(self.sigmas) != list:
self.sigmas = [self.sigmas]

if type(self.num_params) != list:
self.num_params = [self.num_params]

if len(self.sigmas) * len(self.num_params) > 1:
# Inner LCV loop
for num_param in self.num_params:
Expand All @@ -76,14 +76,15 @@ def fit(self, X_train, X_test, alpha_0=None):
for k in range(1,cv+1):
if self.verbose > 0:
print('Training: sigma: %s R: %s' % (sigma, k))
X_test_fold = X_test_shuffled[(k-1)*chunk:k*chunk,:]
j_scores[(num_param,sigma)][k-1] = self._fit(X_train=X_train,
X_test_fold = X_test_shuffled[(k-1)*chunk:k*chunk,:]
j_scores[(num_param,sigma)][k-1] = self._fit(X_train=X_train,
X_test=X_test_fold,
num_parameters = num_param,
sigma=sigma)
j_scores[(num_param,sigma)] = np.mean(j_scores[(num_param,sigma)])

sorted_scores = sorted([x for x in j_scores.iteritems() if np.isfinite(x[1])], key=lambda x :x[1], reverse=True)
sorted_scores = sorted([x for x in j_scores.items() if np.isfinite(x[1])], key=lambda x :x[1], reverse=True)
#print(sorted_scores)
if len(sorted_scores) == 0:
warnings.warn('LCV failed to converge for all values of sigma.')
return self
Expand All @@ -97,59 +98,59 @@ def fit(self, X_train, X_test, alpha_0=None):
self._j = self._fit(X_train=X_train, X_test=X_test_shuffled, num_parameters=self._num_parameters, sigma=self._sigma)

return self # Compatibility with sklearn

def _fit(self, X_train, X_test, num_parameters, sigma, alpha_0=None):
""" Fits the estimator with the given parameters w-hat and returns J"""

num_parameters = num_parameters

if type(num_parameters) == float:
num_parameters = int(X_test.shape[0] * num_parameters)

self._select_param_vectors(X_test=X_test,
self._select_param_vectors(X_test=X_test,
sigma=sigma,
num_parameters=num_parameters)

X_train = self._reshape_X(X_train)
X_test = self._reshape_X(X_test)

if alpha_0 is None:
alpha_0 = np.ones(shape=(num_parameters,1))/float(num_parameters)

self._find_alpha(X_train=X_train,
X_test=X_test,
num_parameters=num_parameters,
epsilon=self.epsilon,
alpha_0 = alpha_0,
sigma=sigma)

return self._calculate_j(X_test,sigma=sigma)

def _calculate_j(self, X_test, sigma):
return np.log(self.predict(X_test,sigma=sigma)).sum()/X_test.shape[0]

def score(self, X_test):
""" Return the J score, similar to sklearn's API """
return self._calculate_j(X_test=X_test, sigma=self._sigma)

@staticmethod
@staticmethod
def _reshape_X(X):
""" Reshape input from mxn to mx1xn to take advantage of numpy broadcasting. """
if len(X.shape) != 3:
return X.reshape((X.shape[0],1,X.shape[1]))
return np.reshape(np.array(X), (X.shape[0],1,X.shape[1]))
return X

def _select_param_vectors(self, X_test, sigma, num_parameters):
""" X_test is the test set. b is the number of parameters. """
""" X_test is the test set. b is the number of parameters. """
indices = np.random.choice(X_test.shape[0], size=num_parameters, replace=False)
self._test_vectors = X_test[indices,:].copy()
self._phi_fitted = True

def _phi(self, X, sigma=None):

if sigma is None:
sigma = self._sigma

if self._phi_fitted:
return np.exp(-np.sum((X-self._test_vectors)**2, axis=-1)/(2*sigma**2))
raise Exception('Phi not fitted.')
Expand All @@ -159,25 +160,23 @@ def _find_alpha(self, alpha_0, X_train, X_test, num_parameters, sigma, epsilon):
b = np.zeros(shape=(num_parameters,1))

A = self._phi(X_test, sigma)
b = self._phi(X_train, sigma).sum(axis=0) / X_train.shape[0]
b = self._phi(X_train, sigma).sum(axis=0) / X_train.shape[0]
b = b.reshape((num_parameters, 1))

out = alpha_0.copy()
for k in range(self.max_iter):
out += epsilon*np.dot(np.transpose(A),1./np.dot(A,out))
out += b*(((1-np.dot(np.transpose(b),out))/np.dot(np.transpose(b),b)))
out = np.maximum(0,out)
out /= (np.dot(np.transpose(b),out))

self._alpha = out
self._fitted = True

def predict(self, X, sigma=None):
""" Equivalent of w(X) from the original paper."""

X = self._reshape_X(X)
if not self._fitted:
raise Exception('Not fitted!')
return np.dot(self._phi(X, sigma=sigma), self._alpha).reshape((X.shape[0],))