Linear regression in R

linear_reg.R

# Kuba 'lambdaofgod' Bartczuk
# linear regression in R 
# implemented with normal equation (check wikipedia)
# 17.08.2014
# TO DO:
# - add normalization
# - add visualization

# loading data
# nothing interesting to see here...
print("Please enter number of features and training examples")
params = scan(file=stdin(), what=integer(), nline=1)
rows = params[2]
cols = params[1] 
X = matrix(0,rows,cols)
y = matrix(0,rows,1)

for(i in 1:rows)
{
    tmp = scan(file = stdin(), nline=1)
    
    X[i,] = tmp[1:cols]
    #normalization should go here
    y[i] = tmp[cols+1]
}


# actual regression
# normal equation:
# theta = (X' * X ) * X' * y
theta = solve_reg(X,y)
print("Theta is: ")
print(theta)

print("Please enter number of test cases")
params = scan(file=stdin(),what=integer(), nline=1)

for (i in 1:params)
    {
        to_test = scan(file = stdin(), nline=1)
        print(to_test %*% theta)
    }


# FUNCTIONS FOR NORMAL EQUATION
# ACTUALLY (KINDA) INTERESTING STUFF

	# kinda obvious, needed for
	pinv_scalar = function(x){
		if (x == 0){
		    return(x)}
		else{
		    return(1/x)}
	    }

	# pseudo-inverse : using SVD pinv(U*D*V') = V*pD*U'
	# where pD is map(pinv_scalar,D) [reciprocated nonzero singular values]
	    pinv = function(M){

		s = svd(M)
		Dinv = make_matrix(pinv_scalar(s$d))
		U = s$u
		V = s$v

		N = V %*% Dinv %*%t(U)

		return(N)
	    }

	# we need this because damn R doesnt get singular values'
	# as matrix but as vector
	    make_matrix = function(v){
		l = length(v)
		N = matrix(0,l,l)

		for (i in 1:l)
		    N[i,i] = v[i]

		return(N)
	    }

	    solve_reg = function(X,y){
		theta = pinv( t(X) %*% X) %*% X %*% y
		return(theta)
	    }