On Tue, Jan 19, 2010 at 4:29 PM, Charles R Harris <charlesr.har...@gmail.com> wrote: > > > On Tue, Jan 19, 2010 at 1:47 PM, Gael Varoquaux > <gael.varoqu...@normalesup.org> wrote: >> >> On Tue, Jan 19, 2010 at 02:22:30PM -0600, Robert Kern wrote: >> > > y = y - np.dot(confounds.T, linalg.lstsq(confounds.T, y)[0]) >> >> > > with np = numpy and linalg = scipy.linalg where scipy calls ATLAS. >> >> > For clarification, are you trying to find the components of the y >> > vectors that are perpendicular to the space spanned by the 10 >> > orthonormal vectors in confounds? >> >> Yes. Actually, what I am doing is calculating partial correlation between >> x and y conditionally to confounds, with the following code: >> >> def cond_partial_cor(y, x, confounds=[]): >> """ Returns the partial correlation of y and x, conditionning on >> confounds. >> """ >> # First orthogonalise y and x relative to confounds >> if len(confounds): >> y = y - np.dot(confounds.T, linalg.lstsq(confounds.T, y)[0]) >> x = x - np.dot(confounds.T, linalg.lstsq(confounds.T, x)[0]) >> return np.dot(x, y)/sqrt(np.dot(y, y)*np.dot(x, x)) >> >> I am not sure that what I am doing is optimal. >> >> > > Most of the time is spent in linalg.lstsq. The length of the vectors >> > > is >> > > 810, and there are about 10 confounds. >> >> > Exactly what are the shapes? y.shape = (810, N); confounds.shape = (810, >> > 10)? >> >> Sorry, I should have been more precise: >> >> y.shape = (810, ) >> confounds.shape = (10, 810) >> > > Column stack the bunch so that the last column is y, then do a qr > decomposition. The last column of q is the (normalized) orthogonal vector > and its amplitude is the last (bottom right) component of r.
do you have to do qr twice, once with x and once with y in the last column or can this be combined? I was trying to do something similar for partial autocorrelation for timeseries but didn't manage or try anything better than repeated leastsq or a variant. > > Chuck > > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > > _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion