probability distributions - problems with the script
2 views (last 30 days)
Show older comments
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);
0 Comments
Accepted Answer
Balaji
on 12 Sep 2023
Edited: Balaji
on 12 Sep 2023
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
0 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!