probability distributions - problems with the script

2 vues (au cours des 30 derniers jours)
ELISABETTA BILLOTTA
ELISABETTA BILLOTTA le 22 Août 2023
Modifié(e) : Balaji le 12 Sep 2023
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);

Réponse acceptée

Balaji
Balaji le 12 Sep 2023
Modifié(e) : Balaji le 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

Plus de réponses (0)

Catégories

En savoir plus sur Response Computation and Visualization dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by