compute_manova.py

import numpy as np

from scipy.stats import f

# to remove, only used for tests
# from sklearn.preprocessing import scale
from .preprocessing_tools import scale

def custom_manova(Y,X,C=None, return_beta = False, bias = None):
    ### /!\ data needs to be scaled, X needs to be (k,1)
    beta = linear_regression(Y, X, C)
    dot_Y = np.dot(Y.T,Y)

    (sign_num, num) = np.linalg.slogdet(dot_Y - np.dot( np.dot(beta, X.T) , np.dot(X, beta.T)))
    # (sign_num, num) = np.linalg.slogdet((Y - X @ beta.T).T@(Y - X @beta.T))
    (sign_denom, denom) = np.linalg.slogdet(dot_Y)
    lamb = np.exp(sign_num*num-(sign_denom*denom))

    fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb

    if return_beta == False :
        return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1)
    else :
        return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1), beta


def custom_mancova(Y,X,C=None):
    ### /!\ data needs to be scaled, X needs to be (k,1)
    print("!!! DEPRECIATED !!! ")
    Y, beta = linear_regression_test(Y, X, C)
    dot_Y = np.dot(Y.T,Y)

    (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
    (sign_denom, denom) = np.linalg.slogdet(dot_Y)
    lamb = np.exp(sign_num*num-(sign_denom*denom))

    fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb
    return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1)



def linear_regression(Y, X, C=None):
    ones = np.ones((X.shape[0],1))
    if C is not None:
        X_C = np.hstack([ones,C,X])
        beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y))
        beta = beta[-1,:].reshape(-1,1) #-1
    else :
        X = np.hstack([ones,X]) #,ones]) #[C,X]
        beta = np.dot(np.linalg.inv(np.dot(X.T,X)) , np.dot(X.T,Y))
        beta = beta[-1,:].reshape(-1,1) 
    return beta


def linear_regression_test(Y, X, C=None):
    ones = np.ones((X.shape[0],1))
    if C is not None:
        X_C = np.hstack([ones,C,X]) #,ones]) #[C,X]
        beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y))
        
        Y = Y- np.dot(X_C[:,:-1],beta[:-1,:]) #Y - X_C[:,:-1] @ beta[:-1,:] # the @ option is slightly slower 13.5 us vs 23 us

        beta = beta[-1,:].reshape(-1,1) #-1
    else :
        beta = np.dot(Y.T,X)/X.shape[0]

    return Y, beta 



from scipy.stats import t

def ttest(beta, X, Y):
    L_p = []

    for i in range(beta.shape[0]) :
        se_b = np.sqrt(1/(X.shape[0]-2)*np.sum((Y[:,i].reshape(-1,1)-beta[i]*X)**2)/np.sum((X-X.mean())**2))
        L_p += [2*(t.sf(np.abs(beta[i])/se_b,X.shape[0]-2))[0]]

    return L_p