Monte carlo simulation code/ unreasonable results in averaging.

Hello,

I am trying to write the following monte carlo simulation code. But, something is wrong and I really need help to fix it.

What I am trying to do is: for every value of k, I am calculating the rate for each one of the relays (MM) and I am selecting the one with the highest rate. For every value of k, I am caluculating evrything for N realizations and I am averaging the value at the end. However, I get the same value for all the values of k and even the value i get is unreasonable! I am sure that i am doing something wrong in the averaging line but I do not know how to fix it!

Any help is appreciated!

function [br]=  Rate_optimum(k)


PL=3;
N=10^5;

MM=3;
si=1;
PT=10.^(5./10);
PR=10.^(0./10);
sign=10.^(0./10);sigd=10.^(0./10);sige=10.^(-7./10);
lamdap=0.1;
 Rate_opt2 =0;

for kk=1:N
     fprintf('Running %d per %d \n',kk,length(k));
   

w_opt=PR;
 
     for bb = 1 : MM   %relays
       f= exprnd(1./lamdap);   
   

 h=exprnd(1./lamdap);  
    g= exprnd(1./lamdap);  

alpha_opt=(k.*si.^(-1))./(PT.*f+(PR.*h.*(PT.*g+sign))+2.*sige);

            SNR1_opt=PT.*f./(sige+sigd);
     SNR2_opt=(w_opt.*h.*g.*PT)./(w_opt.*h.*sign+sige+sigd);
     SNR_SD_opt=SNR1_opt+SNR2_opt;
     Rate=0.5.*(1-alpha_opt).*log2(1+SNR_SD_opt);
                     
            if (Rate  > Rate_opt2 ) 
                Rate_opt2   = Rate ;               
                ID_max        = bb % This is the selected relay from the set of relays we have
            end
     end 
   
     Rate_final=Rate_opt2;

    
     Rate_opt2
end

  
br=   Rate_final/N ;
  br
end

%%%%
%The call function:
% %Test file: run this
% clear all
% close all
% k      = 0:0.5:4; 
% for kk=1:length(k)       
%  br(kk)  =  Rate_optimum(k(kk) );
% end
% 
%  semilogy(k,br,'o','LineWidth',1.5)
% 
%  grid on
 
 

Answer 1

[David Hill]

Very confusing code.

function [br]=  Rate_optimum(k)
PL=3;
N=10^5;
MM=3;
si=1;
PT=10.^(5./10);
PR=10.^(0./10);
sign=10.^(0./10);sigd=10.^(0./10);sige=10.^(-7./10);
lamdap=0.1;
 Rate_opt2 =0;
for kk=1:N
     fprintf('Running %d per %d \n',kk,length(k));
   
w_opt=PR;
 
     for bb = 1 : MM   %relays
       f= exprnd(1./lamdap);   
   
 h=exprnd(1./lamdap);  
    g= exprnd(1./lamdap);  
alpha_opt=(k.*si.^(-1))./(PT.*f+(PR.*h.*(PT.*g+sign))+2.*sige);%k is an array so alpha_out will also be an array
            SNR1_opt=PT.*f./(sige+sigd);
     SNR2_opt=(w_opt.*h.*g.*PT)./(w_opt.*h.*sign+sige+sigd);%everything else is just constant and can be taken outside the loops
     SNR_SD_opt=SNR1_opt+SNR2_opt;
     Rate=0.5.*(1-alpha_opt).*log2(1+SNR_SD_opt);%Rate is also be an array
                     
            if (Rate  > Rate_opt2 ) %does not make sense comparing an array to a scalar
                Rate_opt2   = Rate ;               
                ID_max        = bb % This is the selected relay from the set of relays we have
            end
     end 
   
     Rate_final=Rate_opt2;
    
     Rate_opt2
end
  
br=   Rate_final/N ;
  br
end

How are you differentiating between the different relays? I believe there should be three values for each of your constants (different values for each relay). You never index with bb; therefore, all the relays have the same constants and the rates will be the same. You should be able to do this without any loops.

[Deema42]

Hello,

Thank you for your reply!

the difference between the relays is in the value of the channels: f,g, and h. These are generated randomely so I beleive tha rate for each relay will be different. This is how I differntaite between relays. Will you please suggest a solution for this?

Thank you in advance.

[David Hill]

I believe this is what you are looking for.

N=10^5;
MM=3;
si=1;
PT=10^(5/10);
PR=10^(0/10);
w_opt=PR;
sign=10^(0/10);sigd=10^(0/10);sige=10^(-7/10);
lamdap=0.1;
Rate_opt2 =0;
f=exprnd(1./lamdap,MM,N);
g=exprnd(1./lamdap,MM,N);
h=exprnd(1./lamdap,MM,N);
SNR1_opt=PT*f/(sige+sigd);%MMxNN matrix
SNR2_opt=(w_opt*h.*g*PT)./(w_opt*h*sign+sige+sigd);%MMxNN matrix
SNR_SD_opt=SNR1_opt+SNR2_opt;%MMxNN matrix
alpha_opt=si^(-1)./(PT*f+(PR*h.*(PT*g+sign))+2*sige);%changed alpha_opt so to be constant (MMxNN matrix)
Rate_final=zeros(size(h));
ID_max=zeros(size(h));
for K=1:length(k)%only need to loop the values of k
   Rate=0.5*(1-k(K)*alpha_opt).*log2(1+SNR_SD_opt);%MMxNN matrix
   [Rate_final(K,:)]=max(Rate);%1xNN array
end
[~,ID_max]=max(Rate);%1xNN array (ID_max does not change with k)
br=sum(Rate_final)/N;%1xNN array