<commit_before>import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt
import nibabel as nib

def load_data(f1, f2):
    """
    Return the image data as an array and the convolved time course

    Parameters
    ----------
    f1,f2 : string
        The name of data to process
    Returns
    -------
    tuple:  
        Contains the image data as an array and the convolved time course
    """
    # Load the image as an image object
    img = nib.load('../../data/sub001/BOLD/' + f1 + '.nii.gz')
    # Load the image data as an array
    # Drop the first 4 3D volumes from the array
    data = img.get_data()[..., 4:]
    # Load the pre-written convolved time course
    convolved = np.loadtxt('../../data/convo/' + f2 + '.txt')[4:]
    return(data, convolved)

def reg_voxels_4d(data, convolved):
    design = np.ones((len(convolved), 2))
    design[:, 1] = convolved
    data_2d = np.reshape(data, (-1, data.shape[-1]))
    betas = npl.pinv(design).dot(data_2d.T) #(2, 147456)
    betas_4d = np.reshape(betas.T, data.shape[:-1] + (-1,)) #(64, 64, 36, 2)

    return (betas, betas_4d)

if __name__ == '__main__':
    from sys import argv
    f1 = argv[1]
    f2 = argv[2]
    get_name = f2.replace('/', '_')
    data, convolved = load_data(f1, f2)
    betas_hat, betas_hat_4d = reg_voxels_4d(data, convolved)
    np.savetxt('../../data/beta/' + get_name + '.txt', betas_hat, newline='\r\n')
    plt.imshow(betas_hat_4d[:, :, 14, 0])
    plt.savefig('../../data/beta/task001_run001_conv005.png')
    plt.show()
<commit_msg>add t-test and p-value<commit_after>import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt
import nibabel as nib
from scipy.stats import t as t_dist

def load_data(f1, f2):
    """
    Return the image data as an array and the convolved time course

    Parameters
    ----------
    f1,f2 : string
        The name of data to process
    Returns
    -------
    tuple:  
        Contains the image data as an array and the convolved time course
    """
    # Load the image as an image object
    img = nib.load('../../data/sub001/BOLD/' + f1 + '.nii.gz')
    # Load the image data as an array
    # Drop the first 4 3D volumes from the array
    data = img.get_data()[..., 4:]
    # Load the pre-written convolved time course
    convolved = np.loadtxt('../../data/convo/' + f2 + '.txt')[4:]
    return(data, convolved)

def reg_voxels_4d(data, convolved):
    design = np.ones((len(convolved), 2))
    design[:, 1] = convolved
    data_2d = np.reshape(data, (-1, data.shape[-1]))
    betas = npl.pinv(design).dot(data_2d.T) #(2, 147456)
    betas_4d = np.reshape(betas.T, data.shape[:-1] + (-1,)) #(64, 64, 36, 2)

    return (design, data_2d, betas, betas_4d)   

def RSE(X,Y, betas_hat): #input 2D X, Y and betas_hat
    Y_hat = X.dot(betas_hat)
    res = Y - Y_hat.T
    RSS = np.sum(res ** 2)
    df = X.shape[0] - npl.matrix_rank(X)
    MRSS = RSS / df

    return (MRSS, df) 

def hypothesis(betas_hat, v_cov, df):
    t_stat = betas_hat[:, 1] / np.sqrt(v_cov.diagonal())
    # Get p value for t value suing cumulative density dunction
    ltp = t_dist.cdf(t_stat, df) #lower tail p
    p = 1 - ltp
    return p

if __name__ == '__main__':
    from sys import argv
    f1 = argv[1] #task001_run001/bold
    f2 = argv[2] #task001_run001_conv005
    get_name = f2.replace('/', '_')
    data, convolved = load_data(f1, f2)
    design_mat, data_2d, betas_hat, betas_hat_4d = reg_voxels_4d(data, convolved)
    np.savetxt('../../data/beta/' + get_name + '.txt', betas_hat, newline='\r\n')
    plt.imshow(betas_hat_4d[:, :, 14, 0])
    plt.savefig('../../data/beta/task001_run001_conv005.png')
    #plt.show()
    #get residual standard error
    s2, df = RSE(design_mat, data_2d, betas_hat)

    #get estimator of variance of betas_hat
    v_cov = s2 * npl.inv(design_mat.T.dot(design_mat))

    #T-test on null hypothesis
    p_value = hypothesis(betas_hat, v_cov, df)
    print(p_value)
    np.savetxt(get_name + '_p-value.txt', p_value, newline='\r\n')
{'fraction_non_alphanumeric': 0.1505456349206349, 'fraction_numerical': 0.035962301587301584, 'mean_word_length': 2.7136279926335174, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7121769', 'n_tokens_mistral': 1636, 'n_tokens_neox': 1525, 'n_words': 411}