Графически представить плотность итераций

Вы можете вычислить y для дискретных точек x, установив случайные значения для c и m. Затем с помощью функции hist можно найти "ненормализованную плотность" значений функции для заданного x. Затем вы можете нормализовать его, чтобы получить реальную плотность значений, чтобы общая площадь под кривой распределения составляла в сумме 1.

Чтобы визуализировать это, вы строите сетку сетки [X, Y] по значениям x и y и ставите значения плотности как Z. Теперь вы можете либо построить прибой, либо его контурный график.

Вот код:

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;

Это не очень красиво, но вы можете изменить параметры и поиграть с разбросом/рисунком, чтобы сделать его немного более привлекательным. Также я предположил, что распределение Гаусса вместо вашего случайного (полностью случайное окрашивание даст вам равномерную плотность). Также этот код может быть оптимизирован для скорости.

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

Дает

Это создает гистограмму значений y для каждой позиции в x, затем вычисляет плотность вероятности для каждого значения y и цветов в точках. Здесь значения y для каждой позиции x нормализуются индивидуально, чтобы раскрасить точки в соответствии с общей плотностью графика, который необходимо перенормировать.