Mathematica中文论坛-非官方

标题: 豆粑粑 matlab 画回归散点图 对应 hist图 [打印本页]

---

作者: meatball1982    时间: 2017-1-12 21:40
标题: 豆粑粑 matlab 画回归散点图 对应 hist图
一直想画这样的图。
似乎 matlab自己给出的scatterhist不容易设置想要的形式，特别是左下子图的大小，显示的方式 。
用ax来定义三个图的位置。
用dscatter函数（这是我下载的，调用别人的）
效果还不错。



主函数

1. clear all
   
2. clc
   
3. clf
   
4. 
   
5. h=fun_mm_reg_hist(reg_tr(:,1),reg_tr(:,2),lab_str);
   
6. 
   
7. t = linspace(-1, 1.2, 2000);
   
8. x = (t.^3)+(0.3.*randn(1, 2000));
   
9. y = (t.^3)+(0.3.*randn(1, 2000))+(rand(1,2000)-0.5)*0.01;
   
10. 
    
11. x=(x-min(x))/(max(x)-min(x));
    
12. y=(y-min(y))/(max(y)-min(y));
    
13. h=fun_mm_reg_hist(x',y',{'re','si'});
    
14. fi_na=['../file_imgs/fig_demo'];
    
15. fun_work_li_035_myfig_out(h,fi_na,3);

复制代码

我的函数

1. function [h]=fun_mm_reg_hist(re,si,lab_str)
   
2. 
   
3. n_grid=60;
   
4. 
   
5. x_lin=linspace(0,1,n_grid);
   
6. [x_grid,y_grid]=meshgrid(x_lin,x_lin);
   
7. values = hist3([re si],{x_lin x_lin});
   
8. values = medfilt2(values,[5,5]);
   
9. values(values==0)=NaN;
   
10. 
    
11. [hi_va_tr,hi_bi_tr]=hist(re,linspace(0,1,40));
    
12. hi_va_tr=hi_va_tr./sum(hi_va_tr);
    
13. [hi_va_si,hi_bi_si]=hist(si,linspace(0,1,40));
    
14. hi_va_si=hi_va_si./sum(hi_va_si);
    
15. 
    
16. h=figure(1);
    
17. set(h, 'Position', [1000, 100, 800, 800]);
    
18. ax=axes('position',[0.15 0.15 0.65 0.65]);
    
19. % [val,h]=contourf(x_grid,y_grid,values','edgecolor','none');
    
20. % shading interp
    
21. % caxis([0,max(values(:))]);
    
22. 
    
23. dscatter(re,si)
    
24. % col_mm=flipud(hot);
    
25. col_mm=jet;
    
26. colormap(col_mm)
    
27. % colorbar
    
28. axis equal
    
29. axis xy
    
30. hold on
    
31. 
    
32. % plot(re,si,'b.');
    
33. plot([0 1],[0 1],'linewidth',2);
    
34. view(0,90);
    
35. set(gca,'fontsize',16,'xtick',[0:0.2:1],'ytick',[0:0.2:1])
    
36. xlabel(lab_str{1});
    
37. ylabel(lab_str{2});
    
38. 
    
39. ax=axes('position',[0.15 0.8 0.65 0.15]);
    
40. h_up=bar(hi_bi_tr,hi_va_tr,'b');
    
41. set(gca,'xtick',[],'ytick',[],'xlim',[0 1]);
    
42. 
    
43. ax=axes('position',[0.8 0.15 0.15 0.65]);
    
44. h_right=barh(hi_bi_si,hi_va_si,'b');
    
45. set(gca,'xtick',[],'ytick',[],'ylim',[0 1]);
    
46. 
    
47. h=gcf;

复制代码

下载的函数 dscatter.m

1. function hAxes = dscatter(X,Y, varargin)
   
2. % DSCATTER creates a scatter plot coloured by density.
   
3. %
   
4. %   DSCATTER(X,Y) creates a scatterplot of X and Y at the locations
   
5. %   specified by the vectors X and Y (which must be the same size), colored
   
6. %   by the density of the points.
   
7. %
   
8. %   DSCATTER(...,'MARKER',M) allows you to set the marker for the
   
9. %   scatter plot. Default is 's', square.
   
10. %
    
11. %   DSCATTER(...,'MSIZE',MS) allows you to set the marker size for the
    
12. %   scatter plot. Default is 10.
    
13. %
    
14. %   DSCATTER(...,'FILLED',false) sets the markers in the scatter plot to be
    
15. %   outline. The default is to use filled markers.
    
16. %
    
17. %   DSCATTER(...,'PLOTTYPE',TYPE) allows you to create other ways of
    
18. %   plotting the scatter data. Options are "surf','mesh' and 'contour'.
    
19. %   These create surf, mesh and contour plots colored by density of the
    
20. %   scatter data.
    
21. %
    
22. %   DSCATTER(...,'BINS',[NX,NY]) allows you to set the number of bins used
    
23. %   for the 2D histogram used to estimate the density. The default is to
    
24. %   use the number of unique values in X and Y up to a maximum of 200.
    
25. %
    
26. %   DSCATTER(...,'SMOOTHING',LAMBDA) allows you to set the smoothing factor
    
27. %   used by the density estimator. The default value is 20 which roughly
    
28. %   means that the smoothing is over 20 bins around a given point.
    
29. %
    
30. %   DSCATTER(...,'LOGY',true) uses a log scale for the yaxis.
    
31. %
    
32. %   Examples:
    
33. %
    
34. %       [data, params] = fcsread('SampleFACS');
    
35. %       dscatter(data(:,1),10.^(data(:,2)/256),'log',1)
    
36. %       % Add contours
    
37. %       hold on
    
38. %       dscatter(data(:,1),10.^(data(:,2)/256),'log',1,'plottype','contour')
    
39. %       hold off
    
40. %       xlabel(params(1).LongName); ylabel(params(2).LongName);
    
41. %      
    
42. %   See also FCSREAD, SCATTER.
    
43. 
    
44. % Copyright 2003-2004 The MathWorks, Inc.
    
45. % $Revision:  [        DISCUZ_CODE_5        ]nbsp;  $Date:  $
    
46. 
    
47. % Reference:
    
48. % Paul H. C. Eilers and Jelle J. Goeman
    
49. % Enhancing scatterplots with smoothed densities
    
50. % Bioinformatics, Mar 2004; 20: 623 - 628.
    
51. 
    
52. lambda = [];
    
53. nbins = [];
    
54. plottype = 'scatter';
    
55. contourFlag = false;
    
56. msize = 10;
    
57. marker = 's';
    
58. logy = false;
    
59. filled = true;
    
60. if nargin > 2
    
61.     if rem(nargin,2) == 1
    
62.         error('Bioinfo:IncorrectNumberOfArguments',...
    
63.             'Incorrect number of arguments to %s.',mfilename);
    
64.     end
    
65.     okargs = {'smoothing','bins','plottype','logy','marker','msize','filled'};
    
66.     for j=1:2:nargin-2
    
67.         pname = varargin{j};
    
68.         pval = varargin{j+1};
    
69.         k = strmatch(lower(pname), okargs); %#ok
    
70.         if isempty(k)
    
71.             error('Bioinfo:UnknownParameterName',...
    
72.                 'Unknown parameter name: %s.',pname);
    
73.         elseif length(k)>1
    
74.             error('Bioinfo:AmbiguousParameterName',...
    
75.                 'Ambiguous parameter name: %s.',pname);
    
76.         else
    
77.             switch(k)
    
78.                 case 1  % smoothing factor
    
79.                     if isnumeric(pval)
    
80.                         lambda = pval;
    
81.                     else
    
82.                         error('Bioinfo:InvalidScoringMatrix','Invalid smoothing parameter.');
    
83.                     end
    
84.                 case 2
    
85.                     if isscalar(pval)
    
86.                         nbins = [ pval pval];
    
87.                     else
    
88.                         nbins = pval;
    
89.                     end
    
90.                 case 3
    
91.                     plottype = pval;
    
92.                 case 4
    
93.                     logy = pval;
    
94.                     Y = log10(Y);
    
95.                 case 5
    
96.                     contourFlag = pval;
    
97.                 case 6
    
98.                     marker = pval;
    
99.                 case 7
    
100.                     msize = pval;
     
101.                 case 8
     
102.                     filled = pval;
     
103.             end
     
104.         end
     
105.     end
     
106. end
     
107. 
     
108. minx = min(X,[],1);
     
109. maxx = max(X,[],1);
     
110. miny = min(Y,[],1);
     
111. maxy = max(Y,[],1);
     
112. 
     
113. if isempty(nbins)
     
114.     nbins = [min(numel(unique(X)),200) ,min(numel(unique(Y)),200) ];
     
115. end
     
116. 
     
117. if isempty(lambda)
     
118.     lambda = 20;
     
119. end
     
120. 
     
121. edges1 = linspace(minx, maxx, nbins(1)+1);
     
122. ctrs1 = edges1(1:end-1) + .5*diff(edges1);
     
123. edges1 = [-Inf edges1(2:end-1) Inf];
     
124. edges2 = linspace(miny, maxy, nbins(2)+1);
     
125. ctrs2 = edges2(1:end-1) + .5*diff(edges2);
     
126. edges2 = [-Inf edges2(2:end-1) Inf];
     
127. 
     
128. [n,p] = size(X);
     
129. bin = zeros(n,2);
     
130. % Reverse the columns to put the first column of X along the horizontal
     
131. % axis, the second along the vertical.
     
132. [dum,bin(:,2)] = histc(X,edges1);
     
133. [dum,bin(:,1)] = histc(Y,edges2);
     
134. H = accumarray(bin,1,nbins([2 1])) ./ n;
     
135. G = smooth1D(H,nbins(2)/lambda);
     
136. F = smooth1D(G',nbins(1)/lambda)';
     
137. % = filter2D(H,lambda);
     
138. 
     
139. if logy
     
140.     ctrs2 = 10.^ctrs2;
     
141.     Y = 10.^Y;
     
142. end
     
143. okTypes = {'surf','mesh','contour','image','scatter'};
     
144. k = strmatch(lower(plottype), okTypes); %#ok
     
145. if isempty(k)
     
146.     error('dscatter:UnknownPlotType',...
     
147.         'Unknown plot type: %s.',plottype);
     
148. elseif length(k)>1
     
149.     error('dscatter:AmbiguousPlotType',...
     
150.         'Ambiguous plot type: %s.',plottype);
     
151. else
     
152.     switch(k)
     
153. 
     
154.         case 1 %'surf'
     
155.             fmax=max(F(:));
     
156.              F(F<fmax/100)=NaN;
     
157.             h = surf(ctrs1,ctrs2,F,'edgealpha',0);
     
158.         case 2 % 'mesh'
     
159.             h = mesh(ctrs1,ctrs2,F);
     
160.         case 3 %'contour'
     
161.             [dummy, h] =contour(ctrs1,ctrs2,F);
     
162.         case 4 %'image'
     
163.             nc = 256;
     
164.             F = F./max(F(:));
     
165.             colormap(repmat(linspace(1,0,nc)',1,3));
     
166.             h =image(ctrs1,ctrs2,floor(nc.*F) + 1);
     
167.         case 5 %'scatter'
     
168.             F = F./max(F(:));
     
169.             ind = sub2ind(size(F),bin(:,1),bin(:,2));
     
170.             col = F(ind);
     
171.             if filled
     
172.                 h = scatter(X,Y,msize,col,marker,'filled');
     
173.             else
     
174.                 h = scatter(X,Y,msize,col,marker);
     
175.             end
     
176.     end
     
177. 
     
178. end
     
179. 
     
180. if logy
     
181.     set(gca,'yscale','log');
     
182. end
     
183. if nargout > 0
     
184.     hAxes = get(h,'parent');
     
185. end
     
186. %%%% This method is quicker for symmetric data.
     
187. % function Z = filter2D(Y,bw)
     
188. % z = -1:(1/bw):1;
     
189. % k = .75 * (1 - z.^2);
     
190. % k = k ./ sum(k);
     
191. % Z = filter2(k'*k,Y);
     
192. 
     
193. function Z = smooth1D(Y,lambda)
     
194. [m,n] = size(Y);
     
195. E = eye(m);
     
196. D1 = diff(E,1);
     
197. D2 = diff(D1,1);
     
198. P = lambda.^2 .* D2'*D2 + 2.*lambda .* D1'*D1;
     
199. Z = (E + P) \ Y;
     

复制代码

---

欢迎光临 Mathematica中文论坛-非官方 (http://ilovemathematica.com/)

Powered by Discuz! X3.2