<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">'''
least-squares linear regression
'''

import numpy as np
import scipy.optimize
import itertools
try:
    from rpy2 import robjects
except:
    pass

def nnls(predictors, output):
    '''
    Perform non-negative least squares regression
    to find a linear combination of the columns in predictors
    that predicts output
    
    predictors should be a (n x m) matrix of inputs, where each row
    represents a data point in m dimensions. output should be an 1D
    vector of length n, such that the output attribute of data point
    predictors[i] is in output[i]
    
    Returns: a dict containing the keys 'beta_hat', 'r2', and 'rbar2', the
    fit coeffients (an array of length m), R-squared, and adjusted R-squared respectivly
    
    '''
    if len(predictors.shape) != 2:
        raise TypeError('Predictors must be 2D')
    n, m = predictors.shape
    if len(output.shape) == 2:
        a, b = output.shape
        if a == 1 or b == 1:
            output = np.reshape(output, max(a,b))
        else:
            raise TypeError('output should be either a 1D array or a n x 1 or 1 x n 2D array')
    elif len(output.shape) &gt; 2:
        raise TypeError('output should be either a 1D array or a n x 1 or 1 x n 2D array')
    if n != len(output):
        raise ValueError('Incompatable dimensions')
    
    beta_hat, l2_residual = scipy.optimize.nnls(predictors, output)
    
    # compute r^2 from 
    # https://www.dlitz.net/articles/2008/correlation/
    hat = np.array(np.matrix(predictors) * np.matrix(beta_hat).T)[:,0]
    xy = hat * output
    xy.sort()
    w_xy = np.sum(xy)
    del xy
    u_x = np.sum(np.sort(hat))
    u_y = np.sum(np.sort(output))
    v_x = np.sum(np.sort(hat**2))
    v_y = np.sum(np.sort(output**2))
    r2_stable = (n * w_xy - u_x*u_y)**2 / ((n*v_x - u_x**2) * (n*v_y - u_y**2))
    
    return_value = {'beta_hat': beta_hat, 'r2': r2_stable}#'r2': r2, 'rbar2': rbar2}
    return return_value    


def ols(predictors, output):
    '''Perform ordinary least squares regression
    to find a linear combination of the columns in predictors
    that predicts output
    
    predictors should be a (n x m) matrix of inputs, where each row
    represents a data point in m dimensions. output should be an 1D
    vector of length n, such that the output attribute of data point
    predictors[i] is in output[i]
    
    Returns: a dict containing the keys 'beta_hat', 'r2', and 'rbar2', the
    fit coeffients (an array of length m), R-squared, and adjusted R-squared respectivly
    '''
    if len(predictors.shape) != 2:
        raise TypeError('Predictors must be 2D')
    n, m = predictors.shape
    if len(output.shape) == 2:
        a, b = output.shape
        if a == 1 or b == 1:
            output = np.reshape(output, max(a,b))
        else:
            raise TypeError('output should be either a 1D array or a n x 1 or 1 x n 2D array')
    elif len(output.shape) &gt; 2:
        raise TypeError('output should be either a 1D array or a n x 1 or 1 x n 2D array')
    if n != len(output):
        raise ValueError('Incompatable dimensions')
    
    beta_hat, residuals2, rank, s = np.linalg.lstsq(predictors, output)
    if rank != predictors.shape[1]:
        raise RuntimeError('Predictors is of defficient rank. You have a repeated column?')
    
    ss_err = np.float(residuals2)
    #print ''
    #print np.sum(np.square(output - np.mean(output)))
    #predicted_output = np.zeros_like(output)
    #for i in range(m):
    #    predicted_output += predictors[:,i] * beta_hat[i]
    #print np.sum(np.square(predicted_output - np.mean(output)))
    
    #print 'cov', np.cov(output, predicted_output)[0,1]
    #print 'std*std', np.std(output) * np.std(predicted_output)
    #r = np.cov(output, predicted_output)[0,1] / (np.std(output) * np.std(predicted_output))
    #r2 = np.square(r)
    ss_tot = np.sum(np.square(output - np.mean(output)))
    r2 = 1.0 - (ss_err) / (ss_tot)
    
    # compute r^2 from 
    # https://www.dlitz.net/articles/2008/correlation/
    hat = np.array(np.matrix(predictors) * np.matrix(beta_hat).T)[:,0]
    xy = hat * output
    xy.sort()
    w_xy = np.sum(xy)
    del xy
    u_x = np.sum(np.sort(hat))
    u_y = np.sum(np.sort(output))
    v_x = np.sum(np.sort(hat**2))
    v_y = np.sum(np.sort(output**2))
    r2_stable = (n * w_xy - u_x*u_y)**2 / ((n*v_x - u_x**2) * (n*v_y - u_y**2))
    
    
    output = {'beta_hat': beta_hat, 'r2': r2, 'r2_stable': r2_stable}#'rbar2': rbar2}
    return output
    
def ols_allcomb(predictors, output):
    return_val = {}
    n,m = predictors.shape
    for i in range(1,m):
        for k in itertools.combinations(range(m), i):
            output[k] = ols(predictors[:,k], output)
    return return_val

def nnls_allcomb(predictors, output):
    return_val = {}
    n,m = predictors.shape
    for i in range(1,m+1):
        for k in itertools.combinations(range(m), i):
            return_val[k] = nnls(predictors[:,k], output)
    return return_val

def ols_R(predictors, output):
    r = robjects.r
    n,m = predictors.shape

    print 'Transfering data to R...'
    r['assign']('pred', r('matrix(nrow=%s, ncol=%s)' % (n, m)))
    for i in range(m):
        r['assign']('pred_i', robjects.FloatVector(predictors[:,i]))
        r('pred[,%d] &lt;- pred_i' % (i+1))        
    r['assign']('output', robjects.FloatVector(output))
    print 'Regressing...'
    
    #r('lmout &lt;- lm(output ~ pred + 0)')
    r('lmout &lt;- lm(output ~ pred)')
    #print 'r^2',; r('print(summary(lmout)$r.squared)')
    #print 'coeff'; r('print(summary(lmout)$coefficients[,1])')
    #r('print(summary(lmout))')
    out = r('summary(lmout)')
    
    r('print(summary(lmout))')
    r('library(nnls)')
    r('a = nnls(pred, output)')
    #r('print(influence(lmout))')
    r('library(relaimpo)')
    r('print(calc.relimp(lmout))')
    r('print(confint(lmout, level=0.95))')
    return None

if __name__ == '__main__':
    length = 1000000
    dtype='float64'
    
    X = np.array(np.random.randn(length,3), dtype=dtype)
    y = 2*X[:,0] + X[:,1] + 0.5*X[:,2] + X[:,0] * X[:,1] + np.array(np.random.randn(length), dtype=dtype)
    #predictors = np.random.rand(3,2)
    
    python_result = ols(X,y)
    #print ols(X * ols(X,y)['beta_hat'], y)
    print 'python coeff: %.7f %.7f %.7f' % tuple(python_result['beta_hat'])
    print 'python r2 %.7f\n' % python_result['r2']
    #print 'nnls', nnls(X, y)
    print 'ols', ols(X, y)
    print 'r results\n', ols_R(X, y)</pre></body></html>