%
disp(' ');
disp(' monte_carlo_random.m  ver 1.2  January 14, 2014 ');
disp(' by Tom Irvine  Email: tom@vibrationdata.com ');
disp(' ');
%
clear a;
clear x;
clear random_array;
%
close all;
%
fig_num=1;
%
disp(' Select distribution  ');
disp('  1=uniform  ');
disp('  2=normal   ');
disp('  3=Rayleigh ');
disp('  4=Weibull  ');
disp('  5=Gumbel   ');
disp('  6=Binomial ');
%
idis=input(' ');
disp(' ');
num=input(' Enter number of values to generate:  ');
disp(' ');
a=zeros(num,1);
%
if(idis==1)  % uniform
%
    lower=input(' Enter lower limit  ');
%        
    upper=input(' Enter upper limit  ');
%        
    delta=upper-lower;
%
    for i=1:num
        a(i)=(lower+delta*rand());
    end
%
end
if(idis==2)  % normal 
%
     mean_a=input(' Enter mean ');
    std_dev=input(' Enter std dev ');
    a=mean_a+std_dev*randn(num,1);
%
end
if(idis==3)  % Rayleigh
%
    std_dev=input(' Enter std dev ');
%    
    for i=1:num
%        
        while(1)
%            
           U=rand();
           if(U>1.0e-20)
              break;
           end
%
        end
        a(i)=std_dev*sqrt(-2*log(U));
    end
%
end
if(idis==4)  % Weibull
%
    scale=input(' Enter scale (lambda > 0):  ');
%
    shape=input('      Enter shape (k > 0):  ');
    q=1./shape;
%
    for i=1:num
        while(1)
            r=rand();
            if(r<=0.9999999)
                 break;
            end
        end  
        X= -log(1.-r);
        X=X^q;
        a(i)=scale*X;
    end
%
end
if(idis==5)  % Gumbel
%
    g=-0.577216;
%
    disp(' Select parameter entry pair: ');
    disp('  1= mean, std dev ');
    disp('  2= mode, scale ');
%
    ipair=input(' ');
%
    if(ipair==1)
%        
        mean_a=input(' Enter the mean:  ');
        std_dev=input(' Enter standard deviation: ');
%
        scale=sqrt(6)*std_dev/pi;
        mode=mean_a+g*scale;
%        
    else
        mode=input(' Enter the mode:  ');
       scale=input(' Enter the scale: ');
    end
%
    for i=1:num
        while(1)
           x=rand();
           if(x>1.0e-20)
               break;
           end
        end
        a(i)=mode-scale*log(log(1./x));
    end
%
end
%
if(idis==6)  % Binomial 
%
     n=input(' Enter the number of independent trials:  ');
     p=input(' Enter the probability (0<p<1):  ');
%     
     f=zeros(n,1);
     cumu=zeros(n,1);
     k=0;
     bc = nchoosek(n,k); %% binomial coefficient
%     
     y=(p^k)*((1.-p)^(n-k));
%   
     f(1)=bc*y;
%     
     for k=1:(n-1)
%
       bc = nchoosek(n,k);
       y=(p^k)*((1.-p)^(n-k));
       f(k+1)=bc*y;
%       
     end  
%
     cumu(1)=0.;
     for k=2:n
         cumu(k)=cumu(k-1)+f(k);
     end    
%
     for k=1:num
        u=rand();
%
        for j=1:(n-1)
           if(u>=cumu(j) && u<=cumu(j+1))
              a(k)=j;
              break;
           end
        end
     end
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
mean_a=mean(a);
std_a=std(a);
%
out1=sprintf('\n   mean =%8.4g  \n',mean_a);
out2=sprintf(' 1-sigma=%8.4g  ',std_a);
out3=sprintf(' 2-sigma=%8.4g  ',2*std_a);
out4=sprintf(' 3-sigma=%8.4g  \n',3*std_a);
%
out5=sprintf(' mean + 1-sigma=%8.4g  ',mean_a+std_a);
out6=sprintf(' mean + 2-sigma=%8.4g  ',mean_a+2*std_a);
out7=sprintf(' mean + 3-sigma=%8.4g  \n',mean_a+3*std_a);
%
min_a=min(a);
max_a=max(a);
range=max_a-min_a;
%
out8=sprintf(' min=%8.4g   max=%8.4g   range=%8.4g  \n',min_a,max_a,range);
%
disp(out1);
disp(out2);
disp(out3);
disp(out4);
disp(out5);
disp(out6);
disp(out7);
disp(out8);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
plot(fig_num);
fig_num=fig_num+1;
plot(a);
ylabel('Amplitude');
xlabel('Index');
grid on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Sort data?  1=yes  2=no ');
isd=input(' ');
if(isd==1)
    a=sort(a);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
disp(' Plot histogram?  1=yes  2=no ');
iph=input(' ');
%
if(iph==1)
    [fig_num]=plot_histogram_random(a,fig_num,idis);
end
%
random_array=a;
%
clear a;
%
disp(' ');
disp(' random number store in:  random_array ');