最小二乘法和主成分分析的比较 matlab 儿子的papa

meatball1982

比较ls和pca，实际上，当点比较少的时候，最小二乘法是个不错的选择，但要观察数据的大体方向，pca是个更好的选择。下面说明。

生成在一定范围内的矩形数据，下图中的蓝色的点(我暂时只生成5000个点。)，将这些点进行旋转(这时是-30度),得到图中的红色点。分别用ls和pca分析得到红色点的主要方向。

1. clear all
   
2. clc
   
3. clf
   
4. %% compare lsfit and  pca
   
5. % outline:
   
6. % 01 gen random points in rectangle
   
7. % 02 rotate and shift rectangle
   
8. % 03 ls fit
   
9. % 04 pca
   
10. % 05 plot result
    
11. 
    
12. %% main
    
13. % parameters ---------------------------------------------
    
14. rng default
    
15. n=5000;          % points number
    
16. the=-pi/6;       % rotate angle
    
17. 
    
18. %% 01 gen random points in rectangle ----------------------
    
19. x=rand(n,1)*4-2;
    
20. y=rand(n,1)*3-1.5;
    
21. 
    
22. %% 02 rotate and shift rectangle --------------------------
    
23. x_new= cos(the)*x+sin(the)*y+3;
    
24. y_new=-sin(the)*x+cos(the)*y+5;
    
25. 
    
26. %% 03 ls fit ----------------------------------------------
    
27. p=polyfit(x_new,y_new,1);
    
28. y_p=p(1)*x_new+p(2);
    
29. 
    
30. %% 04 pca -------------------------------------------------
    
31. dat=[x_new';y_new']';
    
32. [coeff,score,latent] = princomp(dat); % new version of matlab, use pca
    
33. 
    
34. % rotate first PC back
    
35. eye_mat=eye(2);
    
36. eye_new=eye_mat*(coeff)^(-1);
    
37. 
    
38. % gen pca line
    
39. x_cen=mean(x_new);
    
40. y_cen=mean(y_new);
    
41. y_p_pca=eye_new(1,2)/eye_new(1,1)*(x_new-x_cen)+y_cen;
    
42. 
    
43. % first PC angle
    
44. atan(eye_new(1,2)/eye_new(1,1))*180/pi
    
45. 
    
46. % second PC angle
    
47. atan(eye_new(2,2)/eye_new(2,1))*180/pi
    
48. 
    
49. %% 05 plot result------------------------------------------
    
50. hold on
    
51. plot(x,y,'b.')                          % plot ori points
    
52. plot(x_new,y_new,'r.')                  % plot rotate points
    
53. plot(x_new,y_p,'k-','linewidth',2)      % ls fit line
    
54. plot(x_new,y_p_pca,'g-','linewidth',2)  % pca line(first PC)
    
55. 
    
56. grid on
    
57. axis equal
    
58. legend('ori point','rot point','ls fit','PCA','Location','southeast')
    
59. 
    
60. % option: output fig
    
61. % h=gcf;
    
62. % fig_na=['../imgs/fig_01_comp_lsfit_pca'];
    
63. % fun_work_li_035_myfig_out(h,fig_na,3);
    
64. 
    
65. %% logs
    
66. % mod :

复制代码