Source code for spm1d.stats.hotellings


# Copyright (C) 2018  Todd Pataky


import numpy as np
from . import _mvbase, _spm

eps        = np.finfo(float).eps   #smallest float, used to avoid divide-by-zero errors



def _T2_onesample_singlenode(y):  #at a single node:
	y        = np.matrix(y)
	n        = y.shape[0]      #nResponses
	m        = y.mean(axis=0)  #mean vector
	W        = np.cov(y.T) + eps    #covariance
	T2       = n * m * np.linalg.inv(W) * m.T
	return float(T2)

def _T2_twosample_singlenode(yA, yB):  #at a single node:
	JA,JB    = yA.shape[0], yB.shape[0]  #nResponses
	yA,yB    = np.matrix(yA), np.matrix(yB)
	mA,mB    = yA.mean(axis=0), yB.mean(axis=0)  #means
	WA,WB    = np.cov(yA.T), np.cov(yB.T)
	W        = ((JA-1)*WA + (JB-1)*WB) / (JA+JB-2) + eps
	T2       = (JA*JB)/float(JA+JB)  * (mB-mA) * np.linalg.inv(W) * (mB-mA).T
	return float(T2)




[docs]def hotellings(Y, mu=None, roi=None): ''' One-sample Hotelling's T2 test. :Parameters: - *Y* --- (J x Q x I) numpy array - *mu* --- scalar or (Q x I) array (datum) :Returns: - T2 : An **spm1d._spm.SPM_T2** instance ''' if Y.ndim==2: if mu is not None: Y = Y - mu T2 = _T2_onesample_singlenode(Y) J,I = Y.shape m,p = float(J)-1, float(I) v1,v2 = p, m return _spm.SPM0D_T2(T2, (v1, v2)) else: if mu is not None: Y = Y - mu nResponses,nNodes,nVectDim = Y.shape T2 = np.array([_T2_onesample_singlenode(Y[:,i,:]) for i in range(nNodes)]) T2 = T2 if roi is None else np.ma.masked_array(T2, np.logical_not(roi)) R = _mvbase._get_residuals_onesample(Y) W = _mvbase._fwhm(R) rCounts = _mvbase._resel_counts(R, W, roi=roi) m,p = float(nResponses)-1, float(nVectDim) v1,v2 = p, m return _spm.SPM_T2(T2, (v1, v2), W, rCounts, None, None, R, roi=roi)
[docs]def hotellings_paired(YA, YB, roi=None): ''' Paired Hotelling's T2 test. :Parameters: - *YA* --- (J x Q x I) numpy array - *YB* --- (J x Q x I) numpy array :Returns: - T2 : An **spm1d._spm.SPM_T2** instance :Note: - A paired Hotelling's test on (YA,YB) is equivalent to a one-sample Hotelling's test on (YB-YA) ''' return hotellings( YB - YA, roi=roi )
[docs]def hotellings2(YA, YB, equal_var=True, roi=None): ''' Two-sample Hotelling's T2 test. :Parameters: - *YA* --- (J x Q x I) numpy array - *YB* --- (J x Q x I) numpy array :Returns: - T2 : An **spm1d._spm.SPM_T2** instance :Note: - Non-sphericity correction not implemented. Equal variance must be assumed by setting "equal_var=True". ''' if equal_var is not True: raise( UserWarning('Non-sphericity correction not implemented. To continue you must assume equal variance and set "equal_var=True".') ) if YA.ndim==2: T2 = _T2_twosample_singlenode(YA, YB) JA,IA = YA.shape JB,IB = YB.shape # v1,v2 = float(IA), float(JA+JB-IA-1) ###incorrect; these are F df, not T2 df v1,v2 = float(IA), float(JA+JB-2) return _spm.SPM0D_T2(T2, (v1, v2)) else: JA,QA,IA = YA.shape JB,QB,IB = YB.shape T2 = np.array([_T2_twosample_singlenode(YA[:,i,:], YB[:,i,:]) for i in range(QA)]) T2 = T2 if roi is None else np.ma.masked_array(T2, np.logical_not(roi)) R = _mvbase._get_residuals_twosample(YA, YB) W = _mvbase._fwhm(R) rCounts = _mvbase._resel_counts(R, W, roi=roi) # v1,v2 = float(IA), float(JA+JB-IA-1) ###incorrect; these are F df, not T2 df v1,v2 = float(IA), float(JA+JB-2) return _spm.SPM_T2(T2, (v1, v2), W, rCounts, None, None, R, roi=roi)