Scilab 中按密度着色的散点图
Scatter plot colored by density in Scilab
我在几列 (table.dat) 中有大量 table 数字数据，我通过
将其作为矩阵导入到 Scilab 6.0
A=fscanfMat('table.dat');
然后将这个矩阵的两列作为平面中点的x坐标和y坐标。命令
scatter(A(:,1),A(:,2),0,".")
现在生成了一个漂亮的点云，但我想根据平面中数据点的数量密度，即附近点的空间密度，为散点图中的每个点着色。例如，高密度区域的点应为深蓝色，低密度区域的点应为红色，中间所有彩虹色的平滑过渡。
在这个帖子中回答了 Python 的问题：
How can I make a scatter plot colored by density in matplotlib?
但这在 Scilab 中如何实现？
您的问题的解决方案是：
- 正在计算您数据的 kernel density estimate (KDE)，
d；
- 使用
rainbowcolormap(n) 创建具有 n 颜色的颜色图 m；
- 像这样绘制数据：
scatter(x,y,s,d,"fill"); set(gcf(),"color_map",m);，其中 s 是图中标记的大小。
因为我无法使用 stixbox toolbox for Scilab，所以我决定想出一个解决这个问题的方法，所以请准备好一个很长的答案。
纯 Scilab 解决方案
首先，我在 Scilab 宏上实现了 kernel_density()。它的输入是 x，一个 n×p 数据矩阵和 h 带宽。它的作用是计算以每个数据点为中心的 circle/sphere/n-sphere 半径 h 内有多少点。
我在统计这个领域不是很有经验，所以我不得不阅读 KDE。事实证明，我的这个解决方案实际上是一种 KDE 方法，它使用带有 constant and equal weight for the neighbors 的内核（因此我将 h 重命名为 "bandwidth" 而不是 "radius"，并且为什么我在计算中添加了一个 2*h*n 因子）。
此外，由于我缺乏知识，我无法实施一种方法来自动为给定的数据集选择最佳 h，因此您必须通过反复试验来选择它。然而，阅读 Scipy implementation of gaussian_kde(), which I saw in the example you provided in your question, and also using hints from this question, and this reference，我想出了一个方法来将可能的 h 的数量减少到 4（如果你的数据有 2 个维度）。也许真正的统计学家可以在评论中验证它，或者提供更好的方法：
- 计算数据集的协方差矩阵；
- 将其平方根乘以斯科特因子：
n ^ (-1 / (p+4))；
- 为所有
h 画图并选择能提供最佳可视化效果的一个。
仍然可以找到原始 kernel_density 函数 here 并且它在大约 10³ 点时工作正常。如果您要处理的不止于此，请继续阅读。
C 实现
如评论部分所述，Scilab 的实施速度相当慢。为了获得更好的结果，我在 C 中实现了 kdec()，并使用 ilib_for_link() 将其链接到 Scilab 宏。但是，这种方法仍然存在问题（请参阅底部的警告说明）。
要在 Scilab 上使用这个函数，你应该有一个兼容的 C 编译器：
- 如果您使用 UNIX 或类 UNIX 系统，则无需担心。
- 如果使用Windows，请按照
mingw toolbox的说明，在执行kde(). 时加载到Scilab环境中
首先，你必须把kdec.c放在当前的Scilab目录中。
//kdec.c
#include <math.h>
void kdec(double f[], double x[], double *h, int *n, int *p){
/* x[]: (n*p)-by-1 array of data
* *h: bandwitdh
* *n: the number of points
* *p: the number of dimensions
* f[]: the output
*
* the local neighborhood density can be defined as (for constant weight):
* f(x0) = sum_from i_to n of K(||x_i - x_0|| <= h) / 2hn
* where: x0 is the observed point, which can have p-dimensions;
* K(a) = {1 if a == True
* {0 if a == False
*/
int n_ = *n; int p_ = *p; double h_ = *h;
int d, j, k;
double dif, norm;
for(j = 0; j < n_; j++){
f[j] = 0;
for(k = 0; k < n_; k++){
norm = 0;
for(d = 0; d < p_; d++){
dif = x[k + d*n_] - x[j + d*n_];
norm = norm + dif * dif;
}
norm = sqrt(norm);
if (norm <= h_){
f[j] = f[j] + 1;
}
}
f[j] = f[j] / (2 * (h_) * (n_));
}
}
然后，设置 kde.sci 调用 kdec C 函数并包装在新的 Scilab kde 函数中。
//kde.sci
if ~isdef('kde') then
ilib_for_link('kdec','kdec.c',[],"c") //compile and create the new shared library
exec('loader.sce',-1); //load library
end
//create a wrapper function to improve interface with interface 'kdec'
function varargout = kde(x,h)
//x: n-by-p matrix of data, each column is a dimension
//h: bandwitdh
[n, p] = size(x); //n: number of points
//p: number of dimensions
x = x(1:$);
if length(h) ~= 1 then
error("kde(x,h): x should be n-by-p matrx; " +...
"h shoud be scalar, positive, and real");
end
f = call('kdec'...
, x , 2, 'd'...
, abs(h), 3, 'd'...
, n , 4, 'i'...
, p , 5, 'i'...
,'out'...
,[n,1] , 1, 'd' );
varargout = list(f)
endfunction
由于我在统计方面没有任何进步，您仍然需要手动设置h。然而，经过多次测试，二维数据的最佳结果似乎是：
scotts_factor = n ^ (-1 / (p+4))
h = sqrt(abs(cov(A))) .* scotts_factor;
h = h(2);
这是一些测试：
exec('kde.sci',-1);
//create data set
n = 1d4;
p = 2;
A = grand((n/2), 1, "nor", 0, 1);
A = [A, A * 3 + grand((n/2), 1, "nor", 0, 1)];
A = [ A ; [ A(:,1) * 0.8 , A(:,2) * 1.3 + 10 ] ];
//calculating bandwidth
scotts_factor = n ^ (-1 / (p+4))
h = sqrt(abs(cov(A))) .* scotts_factor;
h = h(2);
//calculate density
d = kde(A, h);
[d, idx] = gsort(d); //sorting data to plot higher-density points
idx = idx($:-1:1); //over lower-density ones
d = d($:-1:1); //(reversing densities matrix)
A = A(idx,:); //(reordering data matrix)
//plotting
scf(); clf();
scatter(A(:,1), A(:,2), 10, d, "fill");
m = rainbowcolormap(32); //create the rainbow color map
m = m($:-1:1,:); //reverse it to get hotter colors on higher densities
set(gcf(),'color_map',m); //set the desired color map
输出为：
警告说明
即使用C实现后，它仍然是一个高成本的功能。由于两个嵌套的 for 循环，它是 O(n²)。
我做了一些测量，结果如下：
n (points) | 10^3 | 5*10^3 | 10^4 | 10^5
-------------+---------+--------+--------+---------
t (seconds) | 0.13751 | 1.2772 | 4.4545 | 323.34
100k点运行kde()用了5分多钟。既然你说要评估1M点，我也不推荐这个方案。不过，将其与纯 Scilab 解决方案进行比较：后者仅需 5 秒即可处理 10³ 个点（！）。这已经是一个巨大的改进，但恐怕我的解决方案不会变得更好。也许您应该尝试减少样本数量，或者寻找其他计算工具，例如 R.
我在几列 (table.dat) 中有大量 table 数字数据，我通过
将其作为矩阵导入到 Scilab 6.0A=fscanfMat('table.dat');
然后将这个矩阵的两列作为平面中点的x坐标和y坐标。命令
scatter(A(:,1),A(:,2),0,".")
现在生成了一个漂亮的点云，但我想根据平面中数据点的数量密度，即附近点的空间密度，为散点图中的每个点着色。例如，高密度区域的点应为深蓝色，低密度区域的点应为红色，中间所有彩虹色的平滑过渡。
在这个帖子中回答了 Python 的问题： How can I make a scatter plot colored by density in matplotlib?
但这在 Scilab 中如何实现？
您的问题的解决方案是：
- 正在计算您数据的 kernel density estimate (KDE)，
d； - 使用
rainbowcolormap(n)创建具有n颜色的颜色图m； - 像这样绘制数据：
scatter(x,y,s,d,"fill"); set(gcf(),"color_map",m);，其中s是图中标记的大小。
因为我无法使用 stixbox toolbox for Scilab，所以我决定想出一个解决这个问题的方法，所以请准备好一个很长的答案。
纯 Scilab 解决方案
首先，我在 Scilab 宏上实现了 kernel_density()。它的输入是 x，一个 n×p 数据矩阵和 h 带宽。它的作用是计算以每个数据点为中心的 circle/sphere/n-sphere 半径 h 内有多少点。
我在统计这个领域不是很有经验，所以我不得不阅读 KDE。事实证明，我的这个解决方案实际上是一种 KDE 方法，它使用带有 constant and equal weight for the neighbors 的内核（因此我将 h 重命名为 "bandwidth" 而不是 "radius"，并且为什么我在计算中添加了一个 2*h*n 因子）。
此外，由于我缺乏知识，我无法实施一种方法来自动为给定的数据集选择最佳 h，因此您必须通过反复试验来选择它。然而，阅读 Scipy implementation of gaussian_kde(), which I saw in the example you provided in your question, and also using hints from this question, and this reference，我想出了一个方法来将可能的 h 的数量减少到 4（如果你的数据有 2 个维度）。也许真正的统计学家可以在评论中验证它，或者提供更好的方法：
- 计算数据集的协方差矩阵；
- 将其平方根乘以斯科特因子：
n ^ (-1 / (p+4))； - 为所有
h画图并选择能提供最佳可视化效果的一个。
仍然可以找到原始 kernel_density 函数 here 并且它在大约 10³ 点时工作正常。如果您要处理的不止于此，请继续阅读。
C 实现
如评论部分所述，Scilab 的实施速度相当慢。为了获得更好的结果，我在 C 中实现了 kdec()，并使用 ilib_for_link() 将其链接到 Scilab 宏。但是，这种方法仍然存在问题（请参阅底部的警告说明）。
要在 Scilab 上使用这个函数，你应该有一个兼容的 C 编译器：
- 如果您使用 UNIX 或类 UNIX 系统，则无需担心。
- 如果使用Windows，请按照
mingwtoolbox的说明，在执行kde(). 时加载到Scilab环境中
首先，你必须把kdec.c放在当前的Scilab目录中。
//kdec.c
#include <math.h>
void kdec(double f[], double x[], double *h, int *n, int *p){
/* x[]: (n*p)-by-1 array of data
* *h: bandwitdh
* *n: the number of points
* *p: the number of dimensions
* f[]: the output
*
* the local neighborhood density can be defined as (for constant weight):
* f(x0) = sum_from i_to n of K(||x_i - x_0|| <= h) / 2hn
* where: x0 is the observed point, which can have p-dimensions;
* K(a) = {1 if a == True
* {0 if a == False
*/
int n_ = *n; int p_ = *p; double h_ = *h;
int d, j, k;
double dif, norm;
for(j = 0; j < n_; j++){
f[j] = 0;
for(k = 0; k < n_; k++){
norm = 0;
for(d = 0; d < p_; d++){
dif = x[k + d*n_] - x[j + d*n_];
norm = norm + dif * dif;
}
norm = sqrt(norm);
if (norm <= h_){
f[j] = f[j] + 1;
}
}
f[j] = f[j] / (2 * (h_) * (n_));
}
}
然后，设置 kde.sci 调用 kdec C 函数并包装在新的 Scilab kde 函数中。
//kde.sci
if ~isdef('kde') then
ilib_for_link('kdec','kdec.c',[],"c") //compile and create the new shared library
exec('loader.sce',-1); //load library
end
//create a wrapper function to improve interface with interface 'kdec'
function varargout = kde(x,h)
//x: n-by-p matrix of data, each column is a dimension
//h: bandwitdh
[n, p] = size(x); //n: number of points
//p: number of dimensions
x = x(1:$);
if length(h) ~= 1 then
error("kde(x,h): x should be n-by-p matrx; " +...
"h shoud be scalar, positive, and real");
end
f = call('kdec'...
, x , 2, 'd'...
, abs(h), 3, 'd'...
, n , 4, 'i'...
, p , 5, 'i'...
,'out'...
,[n,1] , 1, 'd' );
varargout = list(f)
endfunction
由于我在统计方面没有任何进步，您仍然需要手动设置h。然而，经过多次测试，二维数据的最佳结果似乎是：
scotts_factor = n ^ (-1 / (p+4))
h = sqrt(abs(cov(A))) .* scotts_factor;
h = h(2);
这是一些测试：
exec('kde.sci',-1);
//create data set
n = 1d4;
p = 2;
A = grand((n/2), 1, "nor", 0, 1);
A = [A, A * 3 + grand((n/2), 1, "nor", 0, 1)];
A = [ A ; [ A(:,1) * 0.8 , A(:,2) * 1.3 + 10 ] ];
//calculating bandwidth
scotts_factor = n ^ (-1 / (p+4))
h = sqrt(abs(cov(A))) .* scotts_factor;
h = h(2);
//calculate density
d = kde(A, h);
[d, idx] = gsort(d); //sorting data to plot higher-density points
idx = idx($:-1:1); //over lower-density ones
d = d($:-1:1); //(reversing densities matrix)
A = A(idx,:); //(reordering data matrix)
//plotting
scf(); clf();
scatter(A(:,1), A(:,2), 10, d, "fill");
m = rainbowcolormap(32); //create the rainbow color map
m = m($:-1:1,:); //reverse it to get hotter colors on higher densities
set(gcf(),'color_map',m); //set the desired color map
输出为：
警告说明
即使用C实现后，它仍然是一个高成本的功能。由于两个嵌套的 for 循环，它是 O(n²)。 我做了一些测量，结果如下：
n (points) | 10^3 | 5*10^3 | 10^4 | 10^5
-------------+---------+--------+--------+---------
t (seconds) | 0.13751 | 1.2772 | 4.4545 | 323.34
100k点运行kde()用了5分多钟。既然你说要评估1M点，我也不推荐这个方案。不过，将其与纯 Scilab 解决方案进行比较：后者仅需 5 秒即可处理 10³ 个点（！）。这已经是一个巨大的改进，但恐怕我的解决方案不会变得更好。也许您应该尝试减少样本数量，或者寻找其他计算工具，例如 R.