Scilab中按密度着色的散点图

Posted

技术标签:

【中文标题】Scilab中按密度着色的散点图【英文标题】:Scatter plot colored by density in Scilab 【发布时间】:2018-04-11 21:34:28 【问题描述】:

我有一个包含几列数字的大型数据表 (table.dat),我将其作为矩阵导入到 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 中实现呢?

【问题讨论】:

我想你想看看 Scilab 的 stixbox 工具箱:atoms.scilab.org/toolboxes/stixbox。该 Python 解决方案中的主要功能是来自 SciPy 的 gaussian_kde,并且此工具箱提供了 ksdensity 用于内核平滑密度估计。不幸的是,当我尝试运行它时,我的 Scilab 崩溃了。也许你会更幸运。 【参考方案1】:

通过以下方式解决您的问题:

    计算数据的kernel density estimate (KDE)d; 使用rainbowcolormap(n) 创建颜色映射mn 颜色; 像这样绘制数据:scatter(x,y,s,d,"fill"); set(gcf(),"color_map",m);,其中s 是图中标记的大小。

由于我无法使用stixbox toolbox for Scilab,因此我决定想出一个解决此问题的方法,因此请准备好一个长答案。

纯 Scilab 溶液

首先,我在 Scilab 宏上实现了kernel_density()。它的输入是x,一个n×p 数据矩阵,和h 带宽。它的作用是计算以每个数据点为中心、半径为h 的圆/球/n 球内有多少点。

我在这个统计领域不是很有经验,所以我不得不阅读 KDE。事实证明,我的这个解决方案实际上是一种 KDE 方法,它使用带有constant and equal weight for the neighbors 的内核(因此我将h 重命名为“带宽”而不仅仅是“半径”,以及为什么我添加了一个2*h*n 因子计算)。

另外,由于我缺乏知识,我无法实现为给定数据集自动选择最佳h 的方法,因此您必须通过反复试验来选择它。但是,阅读Scipy implementation of gaussian_kde()(我在您在问题中提供的示例中看到)以及使用来自this question 和this reference 的提示,我想出了一种方法将可能的@ 数量减少到4 987654347@(如果您的数据有 2 个维度)。也许真正的统计学家可以在 cmets 中验证它,或者提供更好的方法:

    计算数据集的协方差矩阵; 将其平方根乘以斯科特因子:n ^ (-1 / (p+4)); 为所有 h 绘制图并选择具有最佳可视化效果的那个。

原来的kernel_density 函数仍然可以在here 中找到,它在大约 10³ 点上运行良好。如果您处理的不止这些,请继续阅读。

C 实现

如 cmets 部分所述,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 

运行kde() 获得 100k 点需要超过 5 分钟。既然你说你想评估 1M 点,我也不推荐这个解决方案。不过,将其与纯 Scilab 解决方案进行比较:后者仅需要 5 秒才能处理 10³ 点(!)。这已经是一个巨大的进步,但恐怕我的解决方案不会变得更好。也许您应该尝试减少样本数量,或者寻找其他计算工具,例如R。

【讨论】:

感谢您的精彩回答!但不幸的是,当它是一个非常大的数据集(如 10^6 点)时,它会永远加载。知道如何加快速度吗? @kolaka 实际上,我也注意到这是一个非常缓慢的解决方案,但我忘了对此进行说明。也许一种解决方案是将其翻译成 C 或 FORTRAN 并使用 call 将 i “导入”到 Scilab 脚本,但我没有对此进行调查。 @kolaka 我更新了这个答案。也许它现在对你有用。 感谢@luispauloml,但由于某种原因,我无法在我的系统(Windows 7)上使用它。我从 gcc 获得多个错误,目前无法修复。不幸的是,我也没有任何在 Scilab 中使用 C 编译器的经验。 @kolaka 处理这些错误是了解调用 C 和 FORTRAN 函数这一特性的好机会。如果你真的想坚持使用 Scilab,也许你应该向某人询问这些错误,或者甚至在 *** 上提出另一个问题。不管怎样,祝你任务顺利,不管它是什么。

以上是关于Scilab中按密度着色的散点图的主要内容,如果未能解决你的问题,请参考以下文章

100天精通Python(可视化篇)——第82天:matplotlib绘制不同种类炫酷散点图参数说明+代码实战(二维散点图三维散点图散点图矩阵)

Matlab中具有密度的散点图

Excel中如何正确地画XY散点图

带有 ggplot2 的散点图按 r 中的特定日期间隔着色

pandas散点图-plot.scatter

r 散点图散点图R.