matlab上的ksdensity只能实现一维的概率密度函数，如果需要二维及以上的概率密度，又该如何操作呢？
核密度计算公式
其中：K(.)为核函数（非负、积分为1，符合概率密度性质，并且均值为0）。有很多种核函数，uniform,triangular, biweight, triweight, Epanechnikov,normal等。
高斯核计算公式如下：
其中，f(x)为位置x处的核密度计算函数，h为距离衰减阈值（即带宽）;n为与位置x的距离小于或等于h的要素点数;k函数则表示空间权重函数。
多维核密度工具包
kde下载地址：https://www.ics.uci.edu/~ihler/code/kde.html
- 下载文件并解压于任意一个文件夹下。这里放在E盘的tmp文件夹下。
- 将” E:\tmp@kde\mex”设为matlab的当前路径
- 命令行输入
>> run('E:\tmp\@kde\mex\makemex.m');
等运行，会出现错误
- 修改错误
相关参考：https://blog.csdn.net/lvfeiya/article/details/72860695
主要有数据类型无法转换、对重载函数调用不明确、未声明变量、不同类型的数据之间无法比较问题。另外，由于修改失误，可能造成无法解析外部符号的错误。
针对数据类型无法转换、对重载函数调用不明确参考上面链接【相关参考】的说明
不同类型的数据之间无法比较问题，修改相关数据类型即可。笔者遇到的是unsigned int和int类型数据无法比较，修改unsigned int成int。如果在修改过程中将函数参数的unsigned int型修改成int，导致无法解析错误，是因为函数参数变了，另外一个cpp文件中对于该函数的声明还是原来的函数，找不到。解决办法：将修改后的函数复制到声明处覆盖掉原来的函数声明即可。
kde的使用
编译成功后，将” E:\tmp”设为matlab的当前路径。
kde结构体构建方式
p = kde(points,ksize,weight,type); % 创建核密度估计类对象
以构建二维概率密度为例，其中points是2行n列的矩阵，ksize是带宽，weight每个样本的权重，type核密度函数类型
例如：**p = kde(data', [100;0.5], [], 'E');根据data创建核密度对象，两个维度的带宽分别为100和0.5，不设置样本权重，核函数类型为E（Epanechnikov)。
得到的p可理解为已经计算出核密度的一整块，但是尺度未知，通过hist(p,[xsize ysize])将这一整块的x划分成xsize份，y划分成ysize份，感官理解为分辨率，划分越细，通过mesh构建的图越细腻。
可选的核函数类型有：
E
G
L
kde可视化案例
%提取符合条件的样本
th=5000;
ind=[11,5];
data=DMLT20M1(:,ind);
sub=find(data(:,1)<th);
data=data(sub,:);
%变量的名称及索引（表格第几列），方便图形命名
varname={'Tu(h)','Td(h)','Telapsed(s)','Dst(km)','MDst(km)'};
index=[5,8,11,12,13];
%最大值，最小值，方便设置坐标轴
maxX=ceil(max(data(:,2)));
minX=ceil(min(data(:,2)));
maxY=ceil(max(data(:,1)));
minY=ceil(min(data(:,1)));
% 创建核密度估计类对象
p = kde(data', [100;0.5], [], 'E');
%可视化p
nx=1000;ny=1000;
z=hist(p,[ny,nx]);
x=0:maxX/nx:maxX;
y=0:maxY/ny:maxY;
X=repmat(x(:,1:end-1),ny,1);
Y=repmat((y(:,1:end-1)'),1,nx);
mesh(X,Y,z);
%图形相关设置
title(['Probability Density of ',char(varname(find(index==ind(2)))),'-', char(varname(find(index==ind(1)))),' in Mode1']);
XLim = get(gca,'XLim');
set(gca,'XMinorTick','on','YMinorTick','on','xlim',[-0.25 24],'ylim',[-0.25 max(max(Y))],'XTick',0:1:24);
xlabel(char(varname(find(index==ind(2)))));
ylabel(char(varname(find(index==ind(1)))));
axis square
h = colorbar;
h.Label.String = 'P';
colormap(jet);
%view(0,90); %加入这句，三维显示变成二维
box on;
grid on;
hold off;
可视化效果
设置view(0,90);
注释掉view(0,90);
网友评论