How to store outputs from Monte Carlo simulation and then calculate probability distribution function?
    5 vues (au cours des 30 derniers jours)
  
       Afficher commentaires plus anciens
    
Hello,
I am beginner with MATLAB so i might have stupid doubts or lack very basic knowledge because this is my first time using matlab.
I need to run MC simulation (1000 times) and store the values of RL_1 , RL_2, and RL_3 (1000 of each) and then based on these values i need to calculate the probability distribution function for each of them, any help?
Here is my code:
clear, close all
freq=csvread('freq.csv',0,0,'A1..A11'); 
freq1 = freq.';
freqchoose = 31.5;
a = 25;
zero = 0;
distance =1:1:25;
distance=distance';
number_of_runs = 1;
Grids1=readtable('SL_S1.csv'); % Read Source level calculated on excel;
Grids2=readtable('SL_S2.csv'); % Read Source level calculated on excel;
Grids3=readtable('SL_S3.csv'); % Read Source level calculated on excel;
u1 = unique(Grids1.id);
u2 = unique(Grids2.id);
u3 = unique(Grids3.id);
for n = 1:number_of_runs %Monte Carlo simualtion
%----------------------------------------------------------------------------------------------------
%SL values per grid cell for all seasons 
for k = 1:numel(u1)
    Grids_season1{k} = Grids1{Grids1.id==u1(k),:}; %Extract submatrices season 1 (per grid cell) from SL_PDF matrix
end
for m = 1:numel(u2)
    Grids_season2{m} = Grids2{Grids2.id==u2(m),:}; %Extract submatrices season 2 (per grid cell) from SL_PDF matrix
end
for n  = 1:numel(u3)
    Grids_season3{n} = Grids3{Grids3.id==u3(n),:}; %Extract submatrices season 3 (per grid cell) from SL_PDF matrix
end
%----------------------------------------------------------------------------------------------------
%Random Source Level for season 1 
RN = rand(1,100); %Random Number
for RR = 1:length(RN)
    ii = find(freq1 == freqchoose); %ii variable is the column number of the freq we chose
    if isempty(ii)
        disp('error')
        return
    end
  for jj = 1:length(Grids_season1)
   SL_freqchoose{jj} = Grids_season1{jj}(:,ii);
  [SL_freqchoose_hist{jj},bincenters{jj}] = hist(SL_freqchoose{jj});
  PDF{jj}=SL_freqchoose_hist{jj}/length(SL_freqchoose{jj});
  BCC{jj} = cumsum(PDF{jj}); %Bincentres Cumulative Sum
  if RN(RR)<BCC{jj}
        SSL_index{jj} = 1;
    else
        SSL_index{jj} = find(BCC{jj}<=RN(RR),1,'last');   
    end
    SSL_1{jj} = bincenters{jj}(SSL_index{jj}); %Synthetic Source Level
  end
%----------------------------------------------------------------------------------------------------
%The probability distribution of the distance between whale and ship -
%season 1
for jj = 1:length(Grids_season1)
distance_cell{jj}=distance;
pdf{jj}(1:25,1)=1/length(distance); %uniform distribution distance from 0 to 25km
BBC1{jj} = cumsum(pdf{jj});
if RN(RR)<BBC1{jj}
        dist_index{jj} = 1;
    else
        dist_index{jj} = find(BBC1{jj}<=RN(RR),1,'last');   
    end
    dist{jj} = (dist_index{jj}); %Synthetic distance value
Pd{jj} = (2*dist{jj}(1,1)*(-4*dist{jj}(1,1)/(a^3)+pi/(a^2)+dist{jj}(1,1)^2/(a^4)));
TL_1{jj} = 10*log10(Pd{jj});
RL_1{jj} = SSL_1{jj} - TL_1{jj};
end
end
%---------------------------------------------------------------------------------------------------
%Random Source Level for season 2
RN = rand(1,100); %Random Number
for RR = 1:length(RN)
    ii = find(freq1 == freqchoose); %ii variable is the column number of the freq we chose
    if isempty(ii)
        disp('error')
        return
    end
  for jj = 1:length(Grids_season2)
   SL_freqchoose{jj} = Grids_season2{jj}(:,ii);
  [SL_freqchoose_hist{jj},bincenters{jj}] = hist(SL_freqchoose{jj});
  PDF{jj}=SL_freqchoose_hist{jj}/length(SL_freqchoose{jj});
  BCC{jj} = cumsum(PDF{jj}); %Bincentres Cumulative Sum
  if RN(RR)<BCC{jj}
        SSL_index{jj} = 1;
    else
        SSL_index{jj} = find(BCC{jj}<=RN(RR),1,'last');   
    end
    SSL_2{jj} = bincenters{jj}(SSL_index{jj}); %Synthetic Source Level
  end
%---------------------------------------------------------------------------------------------------
%The probability distribution of the distance between whale and ship -
%season 2
for jj = 1:length(Grids_season2)
distance_cell{jj}=distance;
pdf{jj}(1:25,1)=1/length(distance); %uniform distribution distance from 0 to 25km
BBC1{jj} = cumsum(pdf{jj});
if RN(RR)<BBC1{jj}
        dist_index{jj} = 1;
    else
        dist_index{jj} = find(BBC1{jj}<=RN(RR),1,'last');   
    end
    dist{jj} = (dist_index{jj}); %Synthetic distance value
Pd{jj} = (2*dist{jj}(1,1)*(-4*dist{jj}(1,1)/(a^3)+pi/(a^2)+dist{jj}(1,1)^2/(a^4)));
TL_2{jj} = 10*log10(Pd{jj});
RL_2{jj} = SSL_2{jj} - TL_2{jj};
end
end
%---------------------------------------------------------------------------------------------------
%Random Source Level for season 3
RN = rand(1,100); %Random Number
for RR = 1:length(RN)
    ii = find(freq1 == freqchoose); %ii variable is the column number of the freq we chose
    if isempty(ii)
        disp('error')
        return
    end
  for jj = 1:length(Grids_season3)
   SL_freqchoose{jj} = Grids_season3{jj}(:,ii);
  [SL_freqchoose_hist{jj},bincenters{jj}] = hist(SL_freqchoose{jj});
  PDF{jj}=SL_freqchoose_hist{jj}/length(SL_freqchoose{jj});
  BCC{jj} = cumsum(PDF{jj}); %Bincentres Cumulative Sum
  if RN(RR)<BCC{jj}
        SSL_index{jj} = 1;
    else
        SSL_index{jj} = find(BCC{jj}<=RN(RR),1,'last');   
    end
    SSL_3{jj} = bincenters{jj}(SSL_index{jj}); %Synthetic Source Level
  end
end
%------------------------------------------------------------------------------------------------------
%The probability distribution of the distance between whale and ship -
%season 3
for jj = 1:length(Grids_season3)
distance_cell{jj}=distance;
pdf{jj}(1:25,1)=1/length(distance); %uniform distribution distance from 0 to 25km
BBC1{jj} = cumsum(pdf{jj});
if RN(RR)<BBC1{jj}
        dist_index{jj} = 1;
    else
        dist_index{jj} = find(BBC1{jj}<=RN(RR),1,'last');   
    end
    dist{jj} = (dist_index{jj}); %Synthetic distance value
Pd{jj} = (2*dist{jj}(1,1)*(-4*dist{jj}(1,1)/(a^3)+pi/(a^2)+dist{jj}(1,1)^2/(a^4)));
TL_3{jj} = 10*log10(Pd{jj});
RL_3{jj} = SSL_3{jj} - TL_3{jj};
end
end 
0 commentaires
Réponses (0)
Voir également
Catégories
				En savoir plus sur Industrial Statistics 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!
