probability distributions - problems with the script

I have this script which I have attached. however, the probability distributions that are produced as results are not correct. in particular the probability values of the "Method 2-mixing" distribution are too low.

to be sure that the script is correct, the script must work even if P1 = 1 and P2 and P3 equal 0.

Can anyone help me fix this script?

A=50; %numero max prove
% N = 80;   % Numero di triplette di n1, n2 e n3 da campionare
N = 150;   % Numero di triplette di n1, n2 e n3 da campionare
noss = 12;  % Numero di osservazioni - Tefra trovati a Monticchio
% noss = 3;  % Numero di osservazioni - Tefra trovati a Ohrid
% n_max = 22;     % Numero massimo di eruzioni del Vesuvio da 22Ka a 1944 tra i 5 e i 30km
P1 = 0.2;   % Probabilità di eruzione di taglia 1 - Papale 2022
P2 = 0.7;   % Probabilità di eruzione di taglia 2 - Papale 2022
P3 = 0.1;   % Probabilità di eruzione di taglia 3 - Papale 2022

p1sito = 0.43;  % Probabilità mappe MA_phi4 taglia 1 a Monticchio
p2sito = 0.64;  % Probabilità mappe MA_phi4 taglia 2 a Monticchio
p3sito = 0.65;  % Probabilità mappe MA_phi4 taglia 3 a Monticchio

% Inizializzare i vari vettori
prob_media_osservazioni_mont = zeros(length(A) - noss + 1, 1);
n1_list = zeros(N, 1);
n2_list = zeros(N, 1);
n3_list = zeros(N, 1);

for n = noss:A
    prob_osservazioni = zeros(N, 1);
    for N_count = 1:N
         for i = 1:N
        % Campiono N triplette di n1, n2, n3 da una multinomiale
        totalEruzioni = mnrnd(n, [P1, P2, P3], N);
        % totalEruzioni = mnrnd(n, [P1, P2, P3]);
        % Distribuisci le eruzioni totali tra n1, n2 e n3 rispettando la distribuzione multinomiale
        n1 = totalEruzioni(1);
        n2 = totalEruzioni(2);
        n3 = totalEruzioni(3);

        % Salva i valori nei rispettivi vettori
        n1_list(i) = n1;
        n2_list(i) = n2;
        n3_list(i) = n3;
         end

        % Step 2: Calcolo delle probabilità delle eruzioni di taglia 1, 2 e 3
        p_k1 = binopdf(0:n1, n1, p1sito);
        p_k2 = binopdf(0:n2, n2, p2sito);
        p_k3 = binopdf(0:n3, n3, p3sito);

        % Step 3: Calcolo la probabilità di avere k1, k2 e k3 dato noss
        probabilities_total = zeros(n1+1, n2+1, n3+1);
        for k1 = 0:n1
            for k2 = 0:n2
                for k3 = 0:n3
                    % Calcola la somma di k1, k2 e k3
                    k_sum = k1 + k2 + k3;

                    % Se la somma non è uguale a n_oss, la probabilità è 0
                    if k_sum ~= noss
                        probabilities_total(k1+1, k2+1, k3+1) = 0;
                    else
                    % Se la somma è uguale a n_oss, la probabilità è il prodotto
                        k_sum = noss;
                        % probabilità dalle binomiali per le tre taglie
                        prob_k1 = p_k1(k1+1);
                        prob_k2 = p_k2(k2+1);
                        prob_k3 = p_k3(k3+1);
                        probabilities_total(k1+1, k2+1, k3+1) = prob_k1 * prob_k2 * prob_k3;
                    end
                end
            end
        end
        % Step 4: Calcolo della probabilità delle osservazioni dati n1, n2,
        % n3 come la somma di tutte le possibili combinazioni di k1, k2, k3
        % che mi danno il numeor di noss(in base al target)
        prob_osservazioni(N_count) = sum(probabilities_total(:));
        % disp(prob_osservazioni);
    end
    
    %Step 5: ripetere il procedimento per le N triplette

    % Step 6: Calcolo la probabilità delle osservazioni come la media di
    % queste probabilità
    prob_media_osservazioni_mont(n - noss + 1) = mean(prob_osservazioni);
end

total_prob_mont = sum(prob_media_osservazioni_mont);
prob_media_osservazioni_mont = prob_media_osservazioni_mont / total_prob_mont;

% Grafico della distribuzione della probabilità - Metodo2
n_values_mont = noss:A;
figure();
bar(n_values_mont, prob_media_osservazioni_mont);
xlabel('da noss ad A');
ylabel('Prob');
title('Metodo 2-mont');

%%%%%%%%%%%%%%%%%%
A_oh=50; %numero max prove
N_oh= 80;   % Numero di triplette di n1, n2 e n3 da campionare
% noss = 12;  % Numero di osservazioni - Tefra trovati a Monticchio
noss_oh = 3;  % Numero di osservazioni - Tefra trovati a Ohrid
% n_max = 22;     % Numero massimo di eruzioni del Vesuvio da 22Ka a 1944 tra i 5 e i 30km
P1_oh = 0.2;   % Probabilità di eruzione di taglia 1 - Papale 2022
P2_oh = 0.7;   % Probabilità di eruzione di taglia 2 - Papale 2022
P3_oh = 0.1;   % Probabilità di eruzione di taglia 3 - Papale 2022

p4sito = 0.31;  % Probabilità mappe MA_phi4 taglia 1 a Ohrid
p5sito = 0.46;  % Probabilità mappe MA_phi4 taglia 2 a Ohrid
p6sito = 0.47;  % Probabilità mappe MA_phi4 taglia 3 a Ohrid
% p4sito = 0.23;  % Probabilità mappe MA_agg taglia 1 a Ohrid
% p5sito = 0.26;  % Probabilità mappe MA_agg taglia 2 a Ohrid
% p6sito = 0.3;  % Probabilità mappe MA_agg taglia 3 a Ohrid

% Inizializzare i vari vettori
prob_media_osservazioni_oh = zeros(length(A_oh) - noss_oh + 1, 1);
n4_list = zeros(N_oh, 1);
n5_list = zeros(N_oh, 1);
n6_list = zeros(N_oh, 1);

for n_oh = noss_oh:A
    prob_osservazioni_oh = zeros(N_oh, 1);
    for N_count_oh = 1:N_oh
         for i_oh = 1:N_oh
        
        totalEruzioni_oh = mnrnd(n_oh, [P1_oh, P2_oh, P3_oh]);
        
        n4 = totalEruzioni_oh(1);
        n5 = totalEruzioni_oh(2);
        n6 = totalEruzioni_oh(3);

        
        n4_list(i_oh) = n4;
        n5_list(i_oh) = n5;
        n6_list(i_oh) = n6;
         end

        
        p_k4 = binopdf(0:n4, n4, p4sito);
        p_k5 = binopdf(0:n5, n5, p5sito);
        p_k6 = binopdf(0:n6, n6, p6sito);

        probabilities_total_oh = zeros(n4+1, n5+1, n6+1);

        for k4 = 0:n4
            for k5 = 0:n5
                for k6 = 0:n6
                    % Calcola la somma di k1, k2 e k3
                    k_sum_oh = k4 + k5 + k6;

                    % Se la somma non è uguale a n_oss, la probabilità è 0
                    if k_sum_oh ~= noss_oh
                        probabilities_total_oh(k4+1, k5+1, k6+1) = 0;
                    else
                    % Se la somma è uguale a n_oss, la probabilità è il prodotto
                        k_sum_oh = noss_oh;
                        % probabilità dalle binomiali per le tre taglie
                        prob_k4 = p_k4(k4+1);
                        prob_k5 = p_k5(k5+1);
                        prob_k6 = p_k6(k6+1);
                        probabilities_total_oh(k4+1, k5+1, k6+1) = prob_k4 * prob_k5 * prob_k6;
                    end
                end
            end
        end
        % Step 4: Calcolo della probabilità delle osservazioni dati n1, n2,
        % n3 come la somma di tutte le possibili combinazioni di k1, k2, k3
        % che mi danno il numeor di noss(in base al target)
        prob_osservazioni_oh(N_count_oh) = sum(probabilities_total_oh(:));
        % disp(prob_osservazioni);
    end
    
    %Step 5: ripetere il procedimento per le N triplette

    % Step 6: Calcolo la probabilità delle osservazioni come la media di
    % queste probabilità
    prob_media_osservazioni_oh(n_oh - noss_oh + 1) = mean(prob_osservazioni_oh);
end

total_prob_oh = sum(prob_media_osservazioni_oh);
prob_media_osservazioni_oh = prob_media_osservazioni_oh / total_prob_oh;


% Grafico della distribuzione della probabilità - Metodo2
n_values_oh = noss_oh:A_oh;
figure();
bar(n_values_oh, prob_media_osservazioni_oh);
xlabel('da noss ad A');
ylabel('Prob');
title('Metodo 2-oh');

%%%%% Mixing tra le multinomiali %%%%%
% Utilizziamo solo le parti dei vettori che corrispondono alle stesse n osservazioni
n_values_mix = noss:max(A, A_oh);
prob_media_osservazioni_mont_mix = prob_media_osservazioni_mont(1:length(n_values_mix));
prob_media_osservazioni_oh_mix = prob_media_osservazioni_oh(1:length(n_values_mix));

% Calcolo della probabilità media ponderata usando le due parti dei vettori
p_media_ponderata = (prob_media_osservazioni_mont_mix + prob_media_osservazioni_oh_mix) / 2;

% Grafico della distribuzione della probabilità - Metodo2 Mixing
figure();
bar(n_values_mix, p_media_ponderata);
xlabel('da noss ad A');
ylabel('Prob');
title('Metodo 2-mixing');
xline(22, 'r', 'LineWidth', 2);

Answer 1

[Balaji]

Hi Elisabetta

I understand that you are facing issues in the code that you have implemented for generating probability distribution for number of volcano eruptions in a period of time.

You can change the for loop for the first distribution as follows:

I have changed the line of code:

totalEruzioni = mnrnd(n, [P1, P2, P3], N);

to:

totalEruzioni = mnrnd(n, [P1, P2, P3]);

and removed the innermost for loop.

for n = noss:A
    prob_osservazioni = zeros(N, 1);
    for N_count = 1:N 
        totalEruzioni = mnrnd(n, [P1, P2, P3]);
        %Distribuisci le eruzioni totali tra n1, n2 e n3 rispettando la distribuzione multinomiale     
        n1 = totalEruzioni(1);
        n2 = totalEruzioni(2);
        n3 = totalEruzioni(3);
        % Step 2: Calcolo delle probabilità delle eruzioni di taglia 1, 2 e 3       
        p_k1 = binopdf(0:n1, n1, p1sito);
        p_k2 = binopdf(0:n2, n2, p2sito);        
        p_k3 = binopdf(0:n3, n3, p3sito);
        % Step 3: Calcolo la probabilità di avere k1, k2 e k3 dato noss        
        probabilities_total = zeros(n1+1, n2+1, n3+1);        
        for k1 = 0:n1           
            for k2 = 0:n2                
                for k3 = 0:n3                    
                    % Calcola la somma di k1, k2 e k3                    
                    k_sum = k1 + k2 + k3;                                        
                    % Se la somma non è uguale a n_oss, la probabilità è 0                    
                    if k_sum ~= noss                        
                        probabilities_total(k1+1, k2+1, k3+1) = 0;                       
                    else                        
                        % Se la somma è uguale a n_oss, la probabilità è il prodotto                        
                        k_sum = noss;                        
                        % probabilità dalle binomiali per le tre taglie                        
                        prob_k1 = p_k1(k1+1);                        
                        prob_k2 = p_k2(k2+1);                        
                        prob_k3 = p_k3(k3+1);                       
                        probabilities_total(k1+1, k2+1, k3+1) = prob_k1 * prob_k2 * prob_k3;                       
                    end                    
                end                
            end            
        end        
        % Step 4: Calcolo della probabilità delle osservazioni dati n1, n2,        
        % n3 come la somma di tutte le possibili combinazioni di k1, k2, k3        
        % che mi danno il numeor di noss(in base al target)        
        prob_osservazioni(N_count) = sum(probabilities_total(:));        
        % disp(prob_osservazioni);        
    end
    %Step 5: ripetere il procedimento per le N triplett    
    % Step 6: Calcolo la probabilità delle osservazioni come la media di    
    % queste probabilità    
    prob_media_osservazioni_mont(n - noss + 1) = mean(prob_osservazioni);    
end

Hope this helps!

Thanks

Balaji