单变量线性回归程序实现

程序流程

1. 读取数据

[plain] view plaincopyprint?

1. data = load('ex1data1.txt');  
2. X = data(:, 1); y = data(:, 2);  
3. m = length(y); % number of training examples  

2. 画散点图 plotData(X,y)

[plain] view plaincopyprint?

1. plot(x,y,'rx','MarkerSize',10);  
2. ylabel('Profit in $10,000s');  
3. xlabel('Population of City in 10,000s');  

3. 梯度下降法

计算代价函数J computerCost(X,y,theta)

[plain] view plaincopyprint?

1. X = [ones(m, 1), data(:,1)]; % Add a column of ones to x  
2. theta = zeros(2, 1) % initialize fitting parameters  

[plain] view plaincopyprint?

1. y1 = X * theta;  
2. J = (y - y1)'*(y - y1)/(2*m);  

确定迭代次数 

[plain] view plaincopyprint?

1. iterations = 1500;  
2. alpha = 0.01;  

迭代θ0和θ1

[plain] view plaincopyprint?

1. for iter = 1:num_iters  
2.   
3.     temp1 = theta(1) - alpha * sum(X*theta-y)/m;  
4.     temp2 = theta(2) - alpha * (X*theta-y)'*X(:,2)/m;  
5.     theta = [temp1; temp2];  
6.      
7.     J_history(iter) = computeCost(X, y, theta);  
8.   
9. end  

4. 画出由θ确定的直线、输出预测数据

[plain] view plaincopyprint?

1. plot(X(:,2), X*theta, '-')  
2. legend('Training data', 'Linear regression')  
3.   
4. % Predict values for population sizes of 35,000 and 70,000  
5. predict1 = [1, 3.5] *theta;  
6. fprintf('For population = 35,000, we predict a profit of %f\n',...  
7.     predict1*10000);  
8. predict2 = [1, 7] * theta;  
9. fprintf('For population = 70,000, we predict a profit of %f\n',...  
10.     predict2*10000);  

5. 可视化代价函数J

[plain] view plaincopyprint?

1. % Grid over which we will calculate J  
2. theta0_vals = linspace(-10, 10, 100);  
3. theta1_vals = linspace(-1, 4, 100);  
4.   
5. % initialize J_vals to a matrix of 0's  
6. J_vals = zeros(length(theta0_vals), length(theta1_vals));  
7.   
8. % Fill out J_vals  
9. for i = 1:length(theta0_vals)  
10.     for j = 1:length(theta1_vals)  
11.       t = [theta0_vals(i); theta1_vals(j)];      
12.       J_vals(i,j) = computeCost(X, y, t);  
13.     end  
14. end  
15.   
16.   
17. % Because of the way meshgrids work in the surf command, we need to   
18. % transpose J_vals before calling surf, or else the axes will be flipped  
19. J_vals = J_vals';  
20. % Surface plot  
21. figure;  
22. surf(theta0_vals, theta1_vals, J_vals)  
23. xlabel('\theta_0'); ylabel('\theta_1');  

6. J的等值线图

[plain] view plaincopyprint?

1. contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))  
2. xlabel('\theta_0'); ylabel('\theta_1');  
3. hold on;  
4. plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);