Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

generates random variates from over 870 univariate distributions

akimpe(A,N,NN)
function [P] = akimpe(A,N,NN)                                       
%             COMPUTES CUBIC POLYNOMIAL INTERPOLATION OF A PERIODIC    
%             FUNCTION OF EQUALLY-SPACED DATA USING AKIMA'S METHOD.    
%             FOR UNEQUALLY SPACED DATA, SEE A.C.M.                    
%             MATH J., V.17, 589 (1970).                               
%             A(N) IS INPUT ARRAY.                                     
%             P(N,4) IS THE RETURNED ARRAY OF POLYNOMIAL COEF:         
%             F(X) = P(I,1) + P(I,2)*(X-X(I)) + P(I,3)*(X-X(I))**2 +   
%                    P(I,4)*(X-X(I))**3                                
%             NN (=N OR =N-5) IS THE NO. OF ACTIVE DATA PTS IN ARRAY A.
%             IF NN = N, THEN THE INPUT ARRAY IS ASSUMED PERIODIC.     
%             IF NN = N-5, THEN THE INPUT ARRAY IS ASSUMED NON-PERIODIC
%             AND MUST BE AUGMENTED WITH FIVE ADDITIONAL DATA PTS AS   
%             FOLLOWS. PTS A(NN+1),A(NN+2), AND A(NN+3) MUST BE        
%             EXTRAPOLATED FROM THE END OF ARRAY A AND PTS A(N), A(N-1)
%             MUST BE EXTRAPOLATED FROM THE BEGINNING OF ARRAY A.      
%             DIMENSION A(N),P(N,4)                                             
%      
      NM1 = N-1;                                                         
      NM2 = N-2;                                                         
      for I = 1:NN 
        NMI = N-I;                                                         
        IP1 = I+1 - N*fix(I/N);                                               
        IP2 = IP1 + 1 - N*fix(I/NM1)*NMI;                                     
        IP3 = IP2 + 1 - N*fix(I/NM2)*NMI*fix((NM1-I)/2);                           
        IM1 = I-1 + N*fix(1/I);                                               
        IM2 = IM1 - 1 + N*fix(2/I)*(I-1);                                     
        S21 = A(IM1) - A(IM2);                                            
        S32 = A(I) - A(IM1);                                               
        S43 = A(IP1) - A(I);                                               
     	S54 = A(IP2) - A(IP1);                                             
        SS21 = S32;                                                        
        SS32 = S43;                                                        
        SS43 = S54;                                                        
        SS54 = A(IP3) - A(IP2);                                            
%             COMPUTE (5-PT) SLOPES, T1 & T2, AT PTS X(I) & X(IP1).    
        AS53 = abs(S54-S43);                                               
        AS31 = abs(S32-S21);                                               
        DENOM = AS53 + AS31;                                               
        T1 =     0.5*(S43+S32);                                            
        if(DENOM >  0.0) 
            T1 = (AS53*S32 + AS31*S43) / DENOM;
        end
        ASS53 = abs(SS54-SS43);                                            
        ASS31 = abs(SS32-SS21);                                            
        DENOMS = ASS53 + ASS31;                                            
        T2 =     0.5*(SS43+SS32);                                          
        if(DENOMS >  0.0) 
            T2 = (ASS53*SS32 + ASS31*SS43) / DENOMS;
        end
%             COMPUTE POLYNOMIAL COEFFICIENTS.                         
        P(I,1) = A(I);                                                     
        P(I,2) = T1;                                                       
        P(I,3) = 3.0*S43 - 2.0*T1 - T2;                                    
        P(I,4) = T1 + T2 - 2.0*S43;                                        
      end                                                          
%                                                                      
      return                                                           
                                                             

Contact us