Dalek
Dalek

Reputation: 4318

Solving linear regression equation for a sparse matrix

I'd like to solve a multivariate linear regression equation for vector X with m elements while I have n observations, Y. If I assume the measurements have Gaussian random errors. How can I solve this problem using python? My problem looks like this:

enter image description here

A simple example of W when m=5, is given as follows:

enter image description here

P.S. I would like to consider the effect of errors and precisely I want to measure the standard deviation of errors.

Upvotes: 1

Views: 1030

Answers (1)

piRSquared
piRSquared

Reputation: 294358

You can do it like this

def myreg(W, Y):
    from numpy.linalg import pinv
    m, n = Y.shape
    k = W.shape[1]
    X = pinv(W.T.dot(W)).dot(W.T).dot(Y)
    Y_hat = W.dot(X)
    Residuals = Y_hat - Y
    MSE = np.square(Residuals).sum(axis=0) / (m - 2)
    X_var = (MSE[:, None] * pinv(W.T.dot(W)).flatten()).reshape(n, k, k)
    Tstat = X / np.sqrt(X_var.diagonal(axis1=1, axis2=2)).T
    return X, Tstat

demo

W = np.array([
        [ 1, -1,  0,  0,  1],
        [ 0, -1,  1,  0,  0],
        [ 1,  0,  0,  1,  0],
        [ 0,  1, -1,  0,  1],
        [ 1, -1,  0,  0, -1],
    ])

Y = np.array([2, 4, 1, 5, 3])[:, None]

X, V = myreg(W, Y)

Upvotes: 2

Related Questions