Logo

Source code for scikits.statsmodels.tools.tools

'''
Utility functions models code
'''

import numpy as np
import numpy.lib.recfunctions as nprf
import numpy.linalg as L
from scipy.interpolate import interp1d
from scipy.linalg import svdvals

def _make_dictnames(tmp_arr, offset=0):
    """
    Helper function to create a dictionary mapping a column number
    to the name in tmp_arr.
    """
    col_map = {}
    for i,col_name in enumerate(tmp_arr):
        col_map.update({i+offset : col_name})
    return col_map

def drop_missing(Y,X=None, axis=1):
    """
    Returns views on the arrays Y and X where missing observations are dropped.

    Y : array-like
    X : array-like, optional
    axis : int
        Axis along which to look for missing observations.  Default is 1, ie.,
        observations in rows.

    Returns
    -------
    Y : array
        All Y where the
    X : array

    Notes
    -----
    If either Y or X is 1d, it is reshaped to be 2d.
    """
    Y = np.asarray(Y)
    if Y.ndim == 1:
        Y = Y[:,None]
    if X is not None:
        X = np.array(X)
        if X.ndim == 1:
            X = X[:,None]
        keepidx = np.logical_and(~np.isnan(Y).any(axis),~np.isnan(X).any(axis))
        return Y[keepidx], X[keepidx]
    else:
        keepidx = ~np.isnan(Y).any(axis)
        return Y[keepidx]

#TODO: needs to better preserve dtype and be more flexible
# ie., if you still have a string variable in your array you don't
# want to cast it to float
#TODO: add name validator (ie., bad names for datasets.grunfeld)
[docs]def categorical(data, col=None, dictnames=False, drop=False): ''' Returns a dummy matrix given an array of categorical variables. Parameters ---------- data : array A structured array, recarray, or array. This can be either a 1d vector of the categorical variable or a 2d array with the column specifying the categorical variable specified by the col argument. col : 'string', int, or None If data is a structured array or a recarray, `col` can be a string that is the name of the column that contains the variable. For all arrays `col` can be an int that is the (zero-based) column index number. `col` can only be None for a 1d array. The default is None. dictnames : bool, optional If True, a dictionary mapping the column number to the categorical name is returned. Used to have information about plain arrays. drop : bool Whether or not keep the categorical variable in the returned matrix. Returns -------- dummy_matrix, [dictnames, optional] A matrix of dummy (indicator/binary) float variables for the categorical data. If dictnames is True, then the dictionary is returned as well. Notes ----- This returns a dummy variable for EVERY distinct variable. If a a structured or recarray is provided, the names for the new variable is the old variable name - underscore - category name. So if the a variable 'vote' had answers as 'yes' or 'no' then the returned array would have to new variables-- 'vote_yes' and 'vote_no'. There is currently no name checking. Examples -------- >>> import numpy as np >>> import scikits.statsmodels.api as sm Univariate examples >>> import string >>> string_var = [string.lowercase[0:5], string.lowercase[5:10], \ string.lowercase[10:15], string.lowercase[15:20], \ string.lowercase[20:25]] >>> string_var *= 5 >>> string_var = np.asarray(sorted(string_var)) >>> design = sm.tools.categorical(string_var, drop=True) Or for a numerical categorical variable >>> instr = np.floor(np.arange(10,60, step=2)/10) >>> design = sm.tools.categorical(instr, drop=True) With a structured array >>> num = np.random.randn(25,2) >>> struct_ar = np.zeros((25,1), dtype=[('var1', 'f4'),('var2', 'f4'), \ ('instrument','f4'),('str_instr','a5')]) >>> struct_ar['var1'] = num[:,0][:,None] >>> struct_ar['var2'] = num[:,1][:,None] >>> struct_ar['instrument'] = instr[:,None] >>> struct_ar['str_instr'] = string_var[:,None] >>> design = sm.tools.categorical(struct_ar, col='instrument', drop=True) Or >>> design2 = sm.tools.categorical(struct_ar, col='str_instr', drop=True) ''' #TODO: add a NameValidator function # catch recarrays and structured arrays if data.dtype.names or data.__class__ is np.recarray: if not col and np.squeeze(data).ndim > 1: raise IndexError("col is None and the input array is not 1d") if isinstance(col, int): col = data.dtype.names[col] # if col is None and len(data.dtype.names) == 1: # col = data.dtype.names[0] tmp_arr = np.unique(data[col]) # if the cols are shape (#,) vs (#,1) need to add an axis and flip _swap = True if data[col].ndim == 1: tmp_arr = tmp_arr[:,None] _swap = False tmp_dummy = (tmp_arr==data[col]).astype(float) if _swap: tmp_dummy = np.squeeze(tmp_dummy).swapaxes(1,0) if not tmp_arr.dtype.names: tmp_arr = np.squeeze(tmp_arr).astype('str').tolist() elif tmp_arr.dtype.names: tmp_arr = np.squeeze(tmp_arr.tolist()).astype('str').tolist() # prepend the varname and underscore, if col is numeric attribute lookup # is lost for recarrays... if col is None: try: col = data.dtype.names[0] except: col = 'var' #TODO: the above needs to be made robust because there could be many # var_yes, var_no varaibles for instance. tmp_arr = [col + '_'+ item for item in tmp_arr] #TODO: test this for rec and structured arrays!!! if drop is True: # if len(data.dtype) is 1 then we have a 1 column array # if len(data.dtype) == 1: if len(data.dtype) <= 1: if tmp_dummy.shape[0] < tmp_dummy.shape[1]: tmp_dummy = np.squeeze(tmp_dummy).swapaxes(1,0) dt = zip(tmp_arr, [tmp_dummy.dtype.str]*len(tmp_arr)) # preserve array type return np.array(map(tuple, tmp_dummy.tolist()), dtype=dt).view(type(data)) data=nprf.drop_fields(data, col, usemask=False, asrecarray=type(data) is np.recarray) data=nprf.append_fields(data, tmp_arr, data=tmp_dummy, usemask=False, asrecarray=type(data) is np.recarray) return data # handle ndarrays and catch array-like for an error elif data.__class__ is np.ndarray or not isinstance(data,np.ndarray): if not isinstance(data, np.ndarray): raise NotImplementedError("Array-like objects are not supported") if isinstance(col, int): offset = data.shape[1] # need error catching here? tmp_arr = np.unique(data[:,col]) tmp_dummy = (tmp_arr[:,np.newaxis]==data[:,col]).astype(float) tmp_dummy = tmp_dummy.swapaxes(1,0) if drop is True: offset -= 1 data = np.delete(data, col, axis=1).astype(float) data = np.column_stack((data,tmp_dummy)) if dictnames is True: col_map = _make_dictnames(tmp_arr, offset) return data, col_map return data elif col is None and np.squeeze(data).ndim == 1: tmp_arr = np.unique(data) tmp_dummy = (tmp_arr[:,None]==data).astype(float) tmp_dummy = tmp_dummy.swapaxes(1,0) if drop is True: if dictnames is True: col_map = _make_dictnames(tmp_arr) return tmp_dummy, col_map return tmp_dummy else: data = np.column_stack((data, tmp_dummy)) if dictnames is True: col_map = _make_dictnames(tmp_arr, offset=1) return data, col_map return data else: raise IndexError("The index %s is not understood" % col) #TODO: add an axis argument to this for sysreg
[docs]def add_constant(data, prepend=False): ''' This appends a column of ones to an array if prepend==False. For ndarrays it checks to make sure a constant is not already included. If there is at least one column of ones then the original array is returned. Does not check for a constant if a structured or recarray is given. Parameters ---------- data : array-like `data` is the column-ordered design matrix prepend : bool True and the constant is prepended rather than appended. Returns ------- data : array The original array with a constant (column of ones) as the first or last column. Notes ----- .. WARNING:: The default of prepend will be changed to True in the next release of statsmodels. We recommend to use an explicit prepend in any permanent code. ''' import warnings warnings.warn("The default of `prepend` will be changed to True in the " "next release, use explicit prepend", FutureWarning) if not data.dtype.names: data = np.asarray(data) if np.any(data[0]==1): ind = np.squeeze(np.where(data[0]==1)) if ind.size == 1 and np.all(data[:,ind] == 1): return data elif ind.size > 1: for col in ind: if np.all(data[:,col] == 1): return data data = np.column_stack((data, np.ones((data.shape[0], 1)))) if prepend: return np.roll(data, 1, 1) else: return_rec = data.__class__ is np.recarray if prepend: ones = np.ones((data.shape[0], 1), dtype=[('const', float)]) data = nprf.append_fields(ones, data.dtype.names, [data[i] for i in data.dtype.names], usemask=False, asrecarray=return_rec) else: data = nprf.append_fields(data, 'const', np.ones(data.shape[0]), usemask=False, asrecarray = return_rec) return data
[docs]def isestimable(C, D): """ From an q x p contrast matrix C and an n x p design matrix D, checks if the contrast C is estimable by looking at the rank of vstack([C,D]) and verifying it is the same as the rank of D. """ if C.ndim == 1: C.shape = (C.shape[0], 1) new = np.vstack([C, D]) if rank(new) != rank(D): return False return True
[docs]def recipr(X): """ Return the reciprocal of an array, setting all entries less than or equal to 0 to 0. Therefore, it presumes that X should be positive in general. """ x = np.maximum(np.asarray(X).astype(np.float64), 0) return np.greater(x, 0.) / (x + np.less_equal(x, 0.))
[docs]def recipr0(X): """ Return the reciprocal of an array, setting all entries equal to 0 as 0. It does not assume that X should be positive in general. """ test = np.equal(np.asarray(X), 0) return np.where(test, 0, 1. / X)
[docs]def clean0(matrix): """ Erase columns of zeros: can save some time in pseudoinverse. """ colsum = np.add.reduce(matrix**2, 0) val = [matrix[:,i] for i in np.flatnonzero(colsum)] return np.array(np.transpose(val))
[docs]def rank(X, cond=1.0e-12): """ Return the rank of a matrix X based on its generalized inverse, not the SVD. """ X = np.asarray(X) if len(X.shape) == 2: D = svdvals(X) return int(np.add.reduce(np.greater(D / D.max(), cond).astype(np.int32))) else: return int(not np.alltrue(np.equal(X, 0.)))
[docs]def fullrank(X, r=None): """ Return a matrix whose column span is the same as X. If the rank of X is known it can be specified as r -- no check is made to ensure that this really is the rank of X. """ if r is None: r = rank(X) V, D, U = L.svd(X, full_matrices=0) order = np.argsort(D) order = order[::-1] value = [] for i in range(r): value.append(V[:,order[i]]) return np.asarray(np.transpose(value)).astype(np.float64) #TODO: sort out the next three classes/functions
class StepFunction: """ A basic step function. Values at the ends are handled in the simplest way possible: everything to the left of x[0] is set to ival; everything to the right of x[-1] is set to y[-1]. Parameters ---------- x : array-like y : array-like ival : float ival is the value given to the values to the left of x[0]. Default is 0. sorted : bool Default is False. Examples -------- >>> import numpy as np >>> from scikits.statsmodels.tools import StepFunction >>> >>> x = np.arange(20) >>> y = np.arange(20) >>> f = StepFunction(x, y) >>> >>> print f(3.2) 3.0 >>> print f([[3.2,4.5],[24,-3.1]]) [[ 3. 4.] [ 19. 0.]] """ def __init__(self, x, y, ival=0., sorted=False): _x = np.asarray(x) _y = np.asarray(y) if _x.shape != _y.shape: raise ValueError('in StepFunction: x and y do not have the same \ shape') if len(_x.shape) != 1: raise ValueError('in StepFunction: x and y must be 1-dimensional') self.x = np.hstack([[-np.inf], _x]) self.y = np.hstack([[ival], _y]) if not sorted: asort = np.argsort(self.x) self.x = np.take(self.x, asort, 0) self.y = np.take(self.y, asort, 0) self.n = self.x.shape[0] def __call__(self, time): tind = np.searchsorted(self.x, time) - 1 _shape = tind.shape return self.y[tind]
[docs]def ECDF(values): """ Return the Empirical CDF of an array as a step function. Parameters ---------- values : array-like Returns ------- Empirical CDF as a step function. """ x = np.array(values, copy=True) x.sort() x.shape = np.product(x.shape,axis=0) n = x.shape[0] y = (np.arange(n) + 1.) / n return StepFunction(x, y)
[docs]def monotone_fn_inverter(fn, x, vectorized=True, **keywords): """ Given a monotone function x (no checking is done to verify monotonicity) and a set of x values, return an linearly interpolated approximation to its inverse from its values on x. """ if vectorized: y = fn(x, **keywords) else: y = [] for _x in x: y.append(fn(_x, **keywords)) y = np.array(y) a = np.argsort(y) return interp1d(y[a], x[a])
[docs]def unsqueeze(data, axis, oldshape): """ Unsqueeze a collapsed array >>> from numpy import mean >>> from numpy.random import standard_normal >>> x = standard_normal((3,4,5)) >>> m = mean(x, axis=1) >>> m.shape (3, 5) >>> m = unsqueeze(m, 1, x.shape) >>> m.shape (3, 1, 5) >>> """ newshape = list(oldshape) newshape[axis] = 1 return data.reshape(newshape)
def chain_dot(*arrs): """ Returns the dot product of the given matrices. Parameters ---------- arrs: argument list of ndarray Returns ------- Dot product of all arguments. Example ------- >>> import numpy as np >>> from scikits.statsmodels.tools import chain_dot >>> A = np.arange(1,13).reshape(3,4) >>> B = np.arange(3,15).reshape(4,3) >>> C = np.arange(5,8).reshape(3,1) >>> chain_dot(A,B,C) array([[1820], [4300], [6780]]) """ return reduce(lambda x, y: np.dot(y, x), arrs[::-1])