以图形方式表示迭代密度

Question

大家好，

我在 MATLAB 中工作，并且正在使用 Monte Carlo 技术来拟合模型。基本上，如果我们假设我的模型是一个简单的函数，例如

y=m*x^2+c

而且我的参数 m 和 c 都在 0.5 和 10 之间变化，我可以从这样的参数空间中随机抽取，每次都获得一个新的 y。如果我绘制我对 y 的所有实现，我会得到类似于下图的东西：

有没有办法表示实现的密度？我的意思是，有没有办法（而不是绘制所有的实现）来获得某种等高线图，它位于我的迭代的最小值和它的颜色代表某个区间内的实现数量的最大值之间？

谢谢大家！

Answer 1

这不是很漂亮，但您可以改变参数并使用散点图/绘图，使其更具视觉吸引力。 我还假设了一个高斯分布而不是你的随机分布（完全随机的颜色会给你一个均匀的密度）。此代码还可以针对速度进行优化。

n = 1000;
l = 100;
x = linspace(1, 10, l);
y = repmat(x.^2, n, 1);
c = repmat(normrnd(1, 1, n, 1), 1, l);
m = repmat(normrnd(1, 1, n, 1), 1, l);

y = y.*m + c;
p = plot(y', '.');
figure; hold on;
for i = 1:l
    [N,edges,bin] = histcounts(y(:, i));
    density = N./sum(N);
    c = zeros(n, 3);
    for j = 1:n
        c(j, :) = [1-density(bin(j))/max(density), 1-density(bin(j))/max(density), 1-density(bin(j))/max(density)];
    end
    scatter(ones(n, 1)*x(i),y(:, i),[],c,'filled'); 
end

给予

这会为 x 中的每个位置创建一个 y 值的直方图，然后计算每个 y 值和点中颜色的概率密度。在这里，每个位置 x 的 y 值都被单独归一化，以根据您需要重新归一化的绘图的整体密度对点进行着色。

Answer 2

您可以为x 的离散点计算y，同时为c 和m 设置随机值。然后使用hist 函数，您可以找到给定x 的函数值的“非标准化密度”。然后，您可以对其进行归一化以获得值的真实密度，以便分布曲线下的总面积总和为1。

为了可视化它，您沿着x 和y 的值构造一个网格[X, Y]，并将密度值设置为Z。现在您可以绘制 surf 或它的等高线图。

代码如下：

clear;

n = 1000000; %number of simulation steps

%parameter ranges
m_min = 0.5; m_max = 10;
c_min = 0.5; c_max = 10;

%x points
x_min = 1; x_max = 4; x_count = 100;
x = linspace(x_min, x_max, x_count);
x2 = x.^2;

y_min = 0; y_max = m_max*x_max*x_max + c_max; y_step = 1;

m = rand(n, 1)*(m_max - m_min) + m_min;
c = rand(n, 1)*(c_max - c_min) + c_min;

c = repmat(c, 1, x_count);

y = m*x2 + c;

x_step = (x_max- x_min)/(x_count-1);
[X, Y] = meshgrid(x_min:x_step:x_max, y_min-y_step:y_step:y_max+y_step);
Z = zeros(size(X));

bins = y_min:y_step:y_max;

for i=1:x_count

    [n_hist,y_hist] = hist(y(:, i), bins);

    %add zeros on both sides to close the profile
    n_hist = [0 n_hist 0];
    y_hist = [y_min-y_step   y_hist   y_max+y_step];

    %normalization
    S = trapz(y_hist,n_hist); %area under the bow
    n_hist = n_hist/S; %scaling of the bow

    Z(:, i) = n_hist';
end

surf(X, Y, Z, 'EdgeColor','none');
colormap jet;
xlim([x_min x_max]);
ylim([y_min y_max]);
xlabel('X');
ylabel('Y');

figure
contour(X,Y,Z, 20);
colormap jet;
colorbar;
grid on;
title('Density as function of X');
xlabel('X');
ylabel('Y');

另一个有趣的视图是每个部分的绘图取决于x 值：

这是这个情节的代码：

clear;

n = 1000000; %number of simulation steps

%parameter ranges
m_min = 0.5; m_max = 10;
c_min = 0.5; c_max = 10;

%x points
x_min = 1; x_max = 4; x_count = 12;
x = linspace(x_min, x_max, x_count);
x2 = x.^2;

m = rand(n, 1)*(m_max - m_min) + m_min;
c = rand(n, 1)*(c_max - c_min) + c_min;

c = repmat(c, 1, x_count);

y = m*x2 + c;

%colors for the plot
colors = ...
[ 0 0 1; 0 1 0; 1 0 0; 0 1 1; 1 0 1; 0 0.75 0.75; 0 0.5 0; 0.75 0.75 0; ...
1 0.50 0.25; 0.75 0 0.75; 0.7 0.7 0.7; 0.8 0.7 0.6; 0.6 0.5 0.4; 1 1 0; 0 0 0 ];

%container for legend entries
legend_list = cell(1, x_count);

for i=1:x_count
    bin_number = 30; %number of histogramm bins
    [n_hist,y_hist] = hist(y(:, i), bin_number);
    n_hist(1) = 0; n_hist(end) = 0; %set first and last values to zero

    %normalization
    S = trapz(y_hist,n_hist); %area under the bow
    n_hist = n_hist/S; %scaling of the bow
    plot(y_hist,n_hist, 'Color', colors(i, :), 'LineWidth', 2);
    hold on;

    legend_list{i} = sprintf('Plot of x = %2.2f', x(i));
end

xlabel('y');
ylabel('pdf(y)');
legend(legend_list);
title('Density depending on x');
grid on;
hold off;