Matlab繪制散點(diǎn)密度圖的教程詳解

更新時(shí)間：2022年02月26日 17:17:03 作者：slandarer

這篇文章主要介紹了如何使用MATLAB繪制散點(diǎn)密度圖(二維核密度)，文中的示例代碼講解詳細(xì)，對(duì)我們學(xué)習(xí)Matlab有一定幫助，需要的可以參考一下

效果

原理也很簡(jiǎn)單，通過matlab自帶的ksdensity獲得網(wǎng)格每一點(diǎn)密度，通過密度擬合曲面，再計(jì)算每個(gè)數(shù)據(jù)點(diǎn)對(duì)應(yīng)的概率，并將概率映射到顏色即可
為了怕大家找不到函數(shù)這次工具函數(shù)放到最前面

1工具函數(shù)完整代碼

function [CData,h,XMesh,YMesh,ZMesh,colorList]=density2C(X,Y,XList,YList,colorList)
[XMesh,YMesh]=meshgrid(XList,YList);
XYi=[XMesh(:) YMesh(:)];
F=ksdensity([X,Y],XYi);
ZMesh=zeros(size(XMesh));
ZMesh(1:length(F))=F;

h=interp2(XMesh,YMesh,ZMesh,X,Y);
if nargin<5
colorList=[0.2700         0    0.3300
    0.2700    0.2300    0.5100
    0.1900    0.4100    0.5600
    0.1200    0.5600    0.5500
    0.2100    0.7200    0.4700
    0.5600    0.8400    0.2700
    0.9900    0.9100    0.1300];
end
colorFunc=colorFuncFactory(colorList);
CData=colorFunc((h-min(h))./(max(h)-min(h)));
colorList=colorFunc(linspace(0,1,100)');

function colorFunc=colorFuncFactory(colorList)
x=(0:size(colorList,1)-1)./(size(colorList,1)-1);
y1=colorList(:,1);y2=colorList(:,2);y3=colorList(:,3);
colorFunc=@(X)[interp1(x,y1,X,'pchip'),interp1(x,y2,X,'pchip'),interp1(x,y3,X,'pchip')];
end
end

2參數(shù)說明

輸入：

X,Y 散點(diǎn)坐標(biāo)
XList,YList 用來構(gòu)造密度曲面網(wǎng)格的序列，其實(shí)就是把XLim，YLim分成小份，例如XList=0:0.1:10
colorList 顏色表mx3數(shù)組，用來構(gòu)造將高度映射到顏色函數(shù)的數(shù)據(jù)表

輸出：

CData各個(gè)點(diǎn)對(duì)應(yīng)顏色
h 各個(gè)點(diǎn)對(duì)應(yīng)核密度
XMesh,YMesh,ZMesh 核密度曲面數(shù)據(jù)
colorList 插值后更細(xì)密的顏色表

3使用方式

假如編寫了如下程序：

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];
scatter(PntSet(:,1),PntSet(:,2),'filled');

結(jié)果：

3.1散點(diǎn)賦色

將上面那段代碼改寫

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];

CData=density2C(PntSet(:,1),PntSet(:,2),-2:0.1:15,-2:0.1:15);
scatter(PntSet(:,1),PntSet(:,2),'filled','CData',CData);

3.2等高線圖

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];


[~,~,XMesh,YMesh,ZMesh,colorList]=density2C(PntSet(:,1),PntSet(:,2),-2:0.1:12,-2:0.1:12);
colormap(colorList)
contourf(XMesh,YMesh,ZMesh,10)

3.3帶直方圖的散點(diǎn)圖

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];

colorList=[0.9400    0.9700    0.9600
    0.8900    0.9300    0.9200
    0.8200    0.9100    0.8800
    0.6900    0.8500    0.7700
    0.5900    0.7800    0.6900
    0.5500    0.7500    0.6500
    0.4500    0.6500    0.5600
    0.4000    0.5800    0.4900
    0.3500    0.5100    0.4200
    0.2500    0.3600    0.3100
    0.1300    0.1700    0.1400];
CData=density2C(PntSet(:,1),PntSet(:,2),-2:0.1:15,-2:0.1:15,colorList);

set(gcf,'Color',[1 1 1]);

% 主分布圖
ax1=axes('Parent',gcf);hold(ax1,'on')
scatter(ax1,PntSet(:,1),PntSet(:,2),'filled','CData',CData);
ax1.Position=[0.1,0.1,0.6,0.6];

% X軸直方圖
ax2=axes('Parent',gcf);hold(ax2,'on')
histogram(ax2,PntSet(:,1),'FaceColor',[0.78 0.88 0.82],...
    'EdgeColor','none','FaceAlpha',0.7)
ax2.Position=[0.1,0.75,0.6,0.15];
ax2.YColor='none';
ax2.XTickLabel='';
ax2.TickDir='out';
ax2.XLim=ax1.XLim;

% Y軸直方圖
ax3=axes('Parent',gcf);hold(ax3,'on')
histogram(ax3,PntSet(:,2),'FaceColor',[0.78 0.88 0.82],...
    'EdgeColor','none','FaceAlpha',0.7,'Orientation','horizontal')
ax3.Position=[0.75,0.1,0.15,0.6];
ax3.XColor='none';
ax3.YTickLabel='';
ax3.TickDir='out';
ax3.YLim=ax1.YLim;

3.4帶直方圖的等高線圖

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];
colorList=[0.9300    0.9500    0.9700
    0.7900    0.8400    0.9100
    0.6500    0.7300    0.8500
    0.5100    0.6200    0.7900
    0.3700    0.5100    0.7300
    0.2700    0.4100    0.6300
    0.2100    0.3200    0.4900
    0.1500    0.2200    0.3500
    0.0900    0.1300    0.2100
    0.0300    0.0400    0.0700];
[~,~,XMesh,YMesh,ZMesh,colorList]=density2C(PntSet(:,1),PntSet(:,2),-2:0.1:13,-2:0.1:13,colorList);


set(gcf,'Color',[1 1 1]);
% 主分布圖
ax1=axes('Parent',gcf);hold(ax1,'on')
colormap(colorList)
contourf(XMesh,YMesh,ZMesh,10,'EdgeColor','none')
ax1.Position=[0.1,0.1,0.6,0.6];
ax1.TickDir='out';

% X軸直方圖
ax2=axes('Parent',gcf);hold(ax2,'on')
[f,xi]=ksdensity(PntSet(:,1));
fill([xi,xi(1)],[f,0],[0.34 0.47 0.71],'FaceAlpha',...
    0.3,'EdgeColor',[0.34 0.47 0.71],'LineWidth',1.2)
ax2.Position=[0.1,0.75,0.6,0.15];
ax2.YColor='none';
ax2.XTickLabel='';
ax2.TickDir='out';
ax2.XLim=ax1.XLim;

% Y軸直方圖
ax3=axes('Parent',gcf);hold(ax3,'on')
[f,yi]=ksdensity(PntSet(:,2));
fill([f,0],[yi,yi(1)],[0.34 0.47 0.71],'FaceAlpha',...
    0.3,'EdgeColor',[0.34 0.47 0.71],'LineWidth',1.2)
ax3.Position=[0.75,0.1,0.15,0.6];
ax3.XColor='none';
ax3.YTickLabel='';
ax3.TickDir='out';
ax3.YLim=ax1.YLim;

4使用方式擴(kuò)展–與ggplot修飾器聯(lián)動(dòng)

ggplot風(fēng)格修飾器：(點(diǎn)擊圖片跳轉(zhuǎn)鏈接)

示例1

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];

ax=gca;
ax.XLim=[-1 13];
ax.YLim=[-1 13];
ax=ggplotAxes2D(ax);

CData=density2C(PntSet(:,1),PntSet(:,2),0:0.1:15,0:0.1:15);
scatter(PntSet(:,1),PntSet(:,2),'filled','CData',CData);

是不是瞬間有那味了：

示例2

PntSet1=mvnrnd([2 3],[1 0;0 2],800);
PntSet2=mvnrnd([6 7],[1 0;0 2],800);
PntSet3=mvnrnd([8 9],[1 0;0 1],800);

PntSet=[PntSet1;PntSet2;PntSet3];

ax=gca;
ax.XLim=[-3 13];
ax.YLim=[-3 13];
ax=ggplotAxes2D(ax);

[~,~,XMesh,YMesh,ZMesh,colorList]=density2C(PntSet(:,1),PntSet(:,2),-2:0.1:12,-2:0.1:12);
colormap(colorList)
contourf(XMesh,YMesh,ZMesh,10)