Python numpy.vander() Examples

def Nordsieck_RKn(self,t0,y,sw0):
        s=self.number_of_steps
        H=(s-1)*self.H
        co_nord=[N.array([1./2,1.]),N.array([2./5,3./5,1.])]
        l=size(y,0)
        y0=y[0,:]
        yf=self.f(t0,y0,sw0)
        
        if l==3:
            co=N.array([co_nord[0]])
            nord_n=N.vander(co_nord[0],self.number_of_steps+1)
            b=y[1:]-y0-co.T*yf
            nord=Sc.solve(nord_n[0:2,0:2],b)
        elif l==4:
            co=N.array([co_nord[1]])
            nord_n=N.vander(co_nord[1],self.number_of_steps+1)
            b=y[1:]-y0-H*co.T*yf
            nord=Sc.solve(nord_n[0:3,0:3],b)
        nord=N.vstack((y0,H*yf,nord[::-1]))       
        return nord 

def linest(*args, **kwargs):  # Excel reference: https://support.office.com/en-us/article/LINEST-function-84d7d0d9-6e50-4101-977a-fa7abf772b6d

    Y = list(args[0].values())
    X = list(args[1].values())

    if len(args) == 3:
        const = args[2]
        if isinstance(const,str):
            const = (const.lower() == "true")
    else:
        const = True

    degree = kwargs.get('degree',1)

    # build the vandermonde matrix
    A = np.vander(X, degree+1)

    if not const:
        # force the intercept to zero
        A[:,-1] = np.zeros((1,len(X)))

    # perform the fit
    (coefs, residuals, rank, sing_vals) = np.linalg.lstsq(A, Y)

    return coefs 

def _ppoly_eval_2(coeffs, breaks, xnew, fill=np.nan):
    """Evaluate piecewise polynomial manually (another way)"""
    a = breaks[0]
    b = breaks[-1]
    K = coeffs.shape[0]

    saveshape = np.shape(xnew)
    xnew = np.ravel(xnew)
    res = np.empty_like(xnew)
    mask = (xnew >= a) & (xnew <= b)
    res[~mask] = fill
    xx = xnew.compress(mask)
    indxs = np.searchsorted(breaks, xx)-1
    indxs = indxs.clip(0, len(breaks))
    pp = coeffs
    diff = xx - breaks.take(indxs)
    V = np.vander(diff, N=K)
    values = np.array([np.dot(V[k, :], pp[:, indxs[k]]) for k in xrange(len(xx))])
    res[mask] = values
    res.shape = saveshape
    return res 

def __call__(self, xnew):
        saveshape = np.shape(xnew)
        xnew = np.ravel(xnew)
        res = np.empty_like(xnew)
        mask = (xnew >= self.a) & (xnew <= self.b)
        res[~mask] = self.fill
        xx = xnew.compress(mask)
        indxs = np.searchsorted(self.breaks, xx)-1
        indxs = indxs.clip(0, len(self.breaks))
        pp = self.coeffs
        diff = xx - self.breaks.take(indxs)
        V = np.vander(diff, N=self.K)
        # values = np.diag(dot(V,pp[:,indxs]))
        values = array([dot(V[k, :], pp[:, indxs[k]]) for k in xrange(len(xx))])
        res[mask] = values
        res.shape = saveshape
        return res 

def reset_ramsey(res, degree=5):
    '''Ramsey's RESET specification test for linear models

    This is a general specification test, for additional non-linear effects
    in a model.


    Notes
    -----
    The test fits an auxiliary OLS regression where the design matrix, exog,
    is augmented by powers 2 to degree of the fitted values. Then it performs
    an F-test whether these additional terms are significant.

    If the p-value of the f-test is below a threshold, e.g. 0.1, then this
    indicates that there might be additional non-linear effects in the model
    and that the linear model is mis-specified.


    References
    ----------
    http://en.wikipedia.org/wiki/Ramsey_RESET_test

    '''
    order = degree + 1
    k_vars = res.model.exog.shape[1]
    #vander without constant and x:
    y_fitted_vander = np.vander(res.fittedvalues, order)[:, :-2] #drop constant
    exog = np.column_stack((res.model.exog, y_fitted_vander))
    res_aux = OLS(res.model.endog, exog).fit()
    #r_matrix = np.eye(degree, exog.shape[1], k_vars)
    r_matrix = np.eye(degree-1, exog.shape[1], k_vars)
    #df1 = degree - 1
    #df2 = exog.shape[0] - degree - res.df_model  (without constant)
    return res_aux.f_test(r_matrix) #, r_matrix, res_aux 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def albrecht_7():
    alpha = 2 * numpy.arange(8) * pi / 8
    s = numpy.array([cos(alpha), sin(alpha)]).T

    alpha = (2 * numpy.arange(8) + 1) * pi / 8
    t = numpy.array([cos(alpha), sin(alpha)]).T

    sqrt21 = sqrt(21)
    wt1, wt2 = (4998 + pm_ * 343 * sqrt21) / 253125
    tau1, tau2 = sqrt((21 - pm_ * sqrt21) / 28)

    # The values are solutions of
    # 4960228*x^4 - 10267740*x^3 + 6746490*x^2 - 1476540*x + 70425 = 0
    sigma2 = roots([4960228, -10267740, 6746490, -1476540, 70425])
    A = numpy.vander(sigma2, increasing=True).T
    b = numpy.array(
        [frac(57719, 675000), frac(9427, 270000), frac(193, 9000), frac(113, 7200)]
    )
    ws = linear_solve(A, b)

    data = [
        (ws[0], sqrt(sigma2[0]) * s),
        (ws[1], sqrt(sigma2[1]) * s),
        (ws[2], sqrt(sigma2[2]) * s),
        (ws[3], sqrt(sigma2[3]) * s),
        (wt1, tau1 * t),
        (wt2, tau2 * t),
    ]

    points, weights = untangle(data)
    return S2Scheme("Albrecht 7", weights, points, 15, _source) 

def albrecht_8():
    alpha = 2 * numpy.arange(10) * pi / 10
    s = numpy.array([cos(alpha), sin(alpha)]).T

    alpha = (2 * numpy.arange(10) + 1) * pi / 10
    t = numpy.array([cos(alpha), sin(alpha)]).T

    m0 = frac(496439663, 13349499975)

    sqrt7 = sqrt(7)
    wt1, wt2 = (125504 + pm_ * 16054 * sqrt7) / 8751645
    tau1, tau2 = sqrt((14 - pm_ * sqrt7) / 18)

    # The values are solutions of
    # 160901628*x^4 - 364759920*x^3 + 274856190*x^2 - 76570340*x
    # + 6054195 = 0
    sigma2 = roots([160901628, -364759920, 274856190, -76570340, 6054195])
    A = numpy.vander(sigma2, increasing=True).T
    b = numpy.array(
        [
            frac(121827491812, 1802182496625),
            frac(48541, 1666980),
            frac(977, 55566),
            frac(671, 52920),
        ]
    )
    ws = linear_solve(A, b)

    data = [
        (m0, z(2)),
        (ws[0], sqrt(sigma2[0]) * s),
        (ws[1], sqrt(sigma2[1]) * s),
        (ws[2], sqrt(sigma2[2]) * s),
        (ws[3], sqrt(sigma2[3]) * s),
        (wt1, tau1 * t),
        (wt2, tau2 * t),
    ]

    points, weights = untangle(data)
    return S2Scheme("Albrecht 8", weights, points, 17, _source) 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.
    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def _setup(self, config):
        torch.manual_seed(config['seed'])
        self.model = ButterflyProduct(size=config['size'],
                                      complex=False,
                                      fixed_order=config['fixed_order'],
                                      softmax_fn=config['softmax_fn'])
        if (not config['fixed_order']) and config['softmax_fn'] == 'softmax':
            self.semantic_loss_weight = config['semantic_loss_weight']
        self.optimizer = optim.Adam(self.model.parameters(), lr=config['lr'])
        self.n_steps_per_epoch = config['n_steps_per_epoch']
        size = config['size']
        # Need to transpose as dct acts on rows of matrix np.eye, not columns
        n = size
        np.random.seed(0)
        x = np.random.randn(n)
        V = np.vander(x, increasing=True)
        self.target_matrix = torch.tensor(V, dtype=torch.float)
        arange_ = np.arange(size)
        dct_perm = np.concatenate((arange_[::2], arange_[::-2]))
        br_perm = bitreversal_permutation(size)
        assert config['perm'] in ['id', 'br', 'dct']
        if config['perm'] == 'id':
            self.perm = torch.arange(size)
        elif config['perm'] == 'br':
            self.perm = br_perm
        elif config['perm'] == 'dct':
            self.perm = torch.arange(size)[dct_perm][br_perm]
        else:
            assert False, 'Wrong perm in config' 

def _setup(self, config):
        torch.manual_seed(config['seed'])
        self.model = ButterflyProduct(size=config['size'],
                                      complex=True,
                                      fixed_order=config['fixed_order'],
                                      softmax_fn=config['softmax_fn'])
        if (not config['fixed_order']) and config['softmax_fn'] == 'softmax':
            self.semantic_loss_weight = config['semantic_loss_weight']
        self.optimizer = optim.Adam(self.model.parameters(), lr=config['lr'])
        self.n_steps_per_epoch = config['n_steps_per_epoch']
        size = config['size']
        n = size
        np.random.seed(0)
        x = np.random.randn(n)
        V = np.vander(x, increasing=True)
        self.target_matrix = torch.tensor(V, dtype=torch.float)
        arange_ = np.arange(size)
        dct_perm = np.concatenate((arange_[::2], arange_[::-2]))
        br_perm = bitreversal_permutation(size)
        assert config['perm'] in ['id', 'br', 'dct']
        if config['perm'] == 'id':
            self.perm = torch.arange(size)
        elif config['perm'] == 'br':
            self.perm = br_perm
        elif config['perm'] == 'dct':
            self.perm = torch.arange(size)[dct_perm][br_perm]
        else:
            assert False, 'Wrong perm in config' 

def matrix_dft(V):
    N = len(V)
    w = np.exp(-2j * np.pi / N)
    col = np.vander([w], N, True)
    W = np.vander(col.flatten(), N, True) / np.sqrt(N)
    return np.dot(W, V) 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def detrend(x, order=1, axis=0):
    """
    Detrend an array with a trend of given order along axis 0 or 1

    Parameters
    ----------
    x : array_like, 1d or 2d
        data, if 2d, then each row or column is independently detrended with the
        same trendorder, but independent trend estimates
    order : int
        specifies the polynomial order of the trend, zero is constant, one is
        linear trend, two is quadratic trend
    axis : int
        axis can be either 0, observations by rows,
        or 1, observations by columns

    Returns
    -------
    detrended data series : ndarray
        The detrended series is the residual of the linear regression of the
        data on the trend of given order.
    """
    if x.ndim == 2 and int(axis) == 1:
        x = x.T
    elif x.ndim > 2:
        raise NotImplementedError('x.ndim > 2 is not implemented until it is needed')

    nobs = x.shape[0]
    if order == 0:
        # Special case demean
        resid = x - x.mean(axis=0)
    else:
        trends = np.vander(np.arange(float(nobs)), N=order + 1)
        beta = np.linalg.pinv(trends).dot(x)
        resid = x - np.dot(trends, beta)

    if x.ndim == 2 and int(axis) == 1:
        resid = resid.T

    return resid 

def albrecht_6():
    # The values are solutions of
    # 11025*x^3 - 19020*x^2 + 9370*x - 1212 = 0
    sigma2 = roots([11025, -19020, 9370, -1212])
    A = numpy.vander(sigma2, increasing=True).T
    b = numpy.array([frac(1432433, 18849024), frac(1075, 31104), frac(521, 25920)])
    B = linear_solve(A, b)

    B0 = frac(2615, 43632)
    C = frac(16807, 933120)

    alpha = 2 * numpy.arange(10) * pi / 10
    rs = numpy.array([cos(alpha), sin(alpha)]).T

    alpha = (2 * numpy.arange(10) + 1) * pi / 10
    uv = numpy.array([cos(alpha), sin(alpha)]).T

    data = [
        (B0, z(2)),
        (B[0], sqrt(sigma2[0]) * rs),
        (B[1], sqrt(sigma2[1]) * rs),
        (B[2], sqrt(sigma2[2]) * rs),
        (C, sqrt(frac(6, 7)) * uv),
    ]

    points, weights = untangle(data)
    return S2Scheme("Albrecht 6", weights, points, 13, _source) 

def compute_expected_coefs(y, degree, with_intercept=True):
    """Helper function to compute expected coefficients from polynomial
    regression"""
    poly_matrix = np.vander(y.index.values, degree + 1)
    if not with_intercept:
        poly_matrix = poly_matrix[:, :-1]
    return np.linalg.lstsq(poly_matrix, y.values, rcond=None)[0] 

def generate_polynomial_series(n, order, coefs=None):
    """Helper function to generate polynomial series of given order and
    coefficients"""
    if coefs is None:
        coefs = np.ones((order + 1, 1))

    x = np.vander(np.arange(n), N=order + 1).dot(coefs)
    return x.ravel() 

def fit_trend(x, order=0):
    """Fit linear regression with polynomial terms of given order

        x : array_like, shape=[n_samples, n_obs]
        Time series data, each sample is fitted separately
    order : int
        The polynomial order of the trend, zero is constant (mean), one is
        linear trend, two is quadratic trend, and so on.

    Returns
    -------
    coefs : ndarray, shape=[n_samples, order + 1]
        Fitted coefficients of polynomial order for each sample, one column
        means order zero, two columns mean order 1
        (linear), three columns mean order 2 (quadratic), etc

    See Also
    -------
    add_trend
    remove_trend
    """
    x = check_array(x)

    if order == 0:
        coefs = np.mean(x, axis=1).reshape(-1, 1)

    else:
        n_obs = x.shape[1]
        index = np.arange(n_obs)
        poly_terms = np.vander(index, N=order + 1)

        # linear least squares fitting using numpy's optimised routine,
        # assuming samples in columns
        # coefs = np.linalg.pinv(poly_terms).dot(x.T).T
        coefs, _, _, _ = np.linalg.lstsq(poly_terms, x.T, rcond=None)

        # returning fitted coefficients in expected format with samples in rows
        coefs = coefs.T

    return coefs 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def test_graph_laplacian():
    for mat in (np.arange(10) * np.arange(10)[:, np.newaxis],
                np.ones((7, 7)),
                np.eye(19),
                np.vander(np.arange(4)) + np.vander(np.arange(4)).T,):
        sp_mat = sparse.csr_matrix(mat)
        for normed in (True, False):
            laplacian = graph_laplacian(mat, normed=normed)
            n_nodes = mat.shape[0]
            if not normed:
                np.testing.assert_array_almost_equal(laplacian.sum(axis=0),
                                                     np.zeros(n_nodes))
            np.testing.assert_array_almost_equal(laplacian.T, laplacian)
            np.testing.assert_array_almost_equal(
                laplacian, graph_laplacian(sp_mat, normed=normed).toarray()) 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def test_graph_laplacian():
    mats = ('np.arange(10) * np.arange(10)[:, np.newaxis]',
            'np.ones((7, 7))',
            'np.eye(19)',
            'sparse.diags([1, 1], [-1, 1], shape=(4,4))',
            'sparse.diags([1, 1], [-1, 1], shape=(4,4)).todense()',
            'np.asarray(sparse.diags([1, 1], [-1, 1], shape=(4,4)).todense())',
            'np.vander(np.arange(4)) + np.vander(np.arange(4)).T',
            )

    for mat_str in mats:
        for normed in (True, False):
            yield _check_graph_laplacian, mat_str, normed 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander 

def vander(x, n=None):
    """
    Masked values in the input array result in rows of zeros.

    """
    _vander = np.vander(x, n)
    m = getmask(x)
    if m is not nomask:
        _vander[m] = 0
    return _vander