function linear_regression(testset, filename, poly, train_method, ss, st)
	%two optional arguments, stepsize (ss) and stopping threshold (st), both real numbers

	%for polynomial fitting
	function A = vandermonde(x, n)
		n = n+1;
		A = ones(size(x, 1), 1);
		for i = 2:n
		  A = [A A(:,i-1).*x];
		end
	end

	%read and prepare data
	X = load(strcat(testset, '/train_x.txt'));
	y = load(strcat(testset, '/train_y.txt'));
	n = size(X,1);
	if poly > 1
		X = vandermonde(X,poly);
	else
		X = [ones(n, 1) X];
	end
	

	%training
	if strcmp('a', train_method)
		%analytic solution
		B = inv(X'*X)*X'*y;
	elseif strcmp('g', train_method)
		%gradient descent solution
		oldloss = realmax;
		B = zeros(size(X,2), 1);
		rr = realmax;
		while rr > st
			B = B - 2/n*ss*X'*(X*B - y);
			newloss = (1/n)*sum((X*B - y).^2);
			rr = (oldloss - newloss)/oldloss;
			oldloss = newloss;
		end
	end

	%predict
	devX = load(strcat(testset, '/dev_x.txt'));
	devy = load(strcat(testset, '/dev_y.txt'));
	if poly > 1
		devX = vandermonde(devX, poly);
	else
		devX = [ones(size(devX, 1), 1) devX];
	end
	predictions = devX*B;
	%evaluate print to standard out
	fprintf('%.3e\n', (1/size(devy,1))*sum((predictions - devy).^2));

	%print results
	f = fopen(strcat(testset,'/',train_method,filename,'.model'), 'w');
	fprintf(f, '%.3e ', B);
	f = fopen(strcat(testset,'/',train_method,filename,'.pred'), 'w');
	fprintf(f, '%.3e\n', predictions');
end
