Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

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