Friday, January 29, 2016

Special functions: Kummer confluent hypergeometric functions



All the elements presented are not warranted to be correct or free from defects. 
Please report any errors found to  afstblogs@gmail.com



For a wide set of formulas for confluent hypergeometric functions access http://dlmf.nist.gov/13. here we will examine only confluent hypergeometric functions of the first and second kind that are the solution of the Kummer differential equation. Wikipedia presents the following formulas:




Kummer's function of the first kind


This function is also presented as 1F1.

Kummer's function of the second kind



The function of the first kind will be used soon as part of the expression of the probability density functions of the noncentral Student's t and sample correlation distributions.



Excel formulas and user defined functions


The only function presented is the one for the firs kind contained in the function XN.xlam6056M add-in that is free and can be downloaded from http://www.bowdoin.edu/~rdelevie/excellaneous/.  

It is not straightforward to compute this function without using VBA. In alternative the computation with Maple Excel add-in is presented. In future messages its use will be explained.


XNumbers60


'ATTENTION 
'Const TOL# = 1.6 * Ten_16
    'changed to
'Const TOL# = 1.6 * 1E-16

'return the M(a,b,z) Kummer confluent hypergeometric function of 1st kind
'it uses the serie expansion method
'
Function Kummer1(a, b, z)
Dim n&, kn#, PN#, dY#, y#, ua#
'Const TOL# = 1.6 * Ten_16
Const TOL# = 1.6 * 1E-16
kn = 1#: PN = 1#:  y = 1#
'Kummer's transformation for negative value of z
If z < 0 Then ua = b - a Else ua = a
For n = 1 To 999
    kn = kn * (ua + n - 1) / (b + n - 1)
    PN = PN * Abs(z) / n
    dY = kn * PN
    y = y + dY
    If Abs(dY) < TOL * Abs(y) Then Exit For
Next n
If n < 999 Then
    If z < 0 Then y = y * Exp(z)  'Kummer's anti-transformation
    Kummer1 = y
Else
    Kummer1 = "?"
End If
End Function

Matlab

Matlab has not a function to compute this function, meaning that a user defied must be used. The one presented below is from
http://www.mathworks.com/matlabcentral/fileexchange/29766-confluent-hypergeometric-function/content/kummer.m


function f = kummer(a,b,x)

% This function estimates the Kummer function with the specified tolerance
% the generalized hypergeometric series, noted below.  This solves Kummer's
% differential equation:
%
%       x*g''(x) + (b - x)*g'(x) - a*g(x) = 0

% Default tolerance is tol = 1e-10.  Feel free to change this as needed.
tol = 1e-10;

% Estimates the value by summing powers of the generalized hypergeometric
% series:
%
%       sum(n=0-->Inf)[(a)_n*x^n/{(b)_n*n!}
%
% until the specified tolerance is acheived.

term = x*a/b;
f = 1 + term;
n = 1;
an = a;
bn = b;
nmin = 10;
while(n < nmin)||max(abs(term) > tol)
  n = n + 1;
  an = an + 1;
  bn = bn + 1;
  term = x.*term*an/bn/n;
  f = f + term;
end

% VERSION INFORMATION
% v1 - Written to support only scalar inputs for x
% v2 - Changed to support column inputs for x by using the repmat
% command and using matrix multiplication to achieve the desired sum
%
% v3 - Credit goes to Ben Petschel for making this suggestion.
%    The previous method of creating vectors for multiplication to
%    produce the sum was replaced by a while loop that executes
%    until a certain tolerance is achieved.  My previous thinking
%    was avoiding a loop would produce a code that would execute
%    faster.  Ben pointed out this is not necessarily true, and not
%    true in this case.  Not only does the while loop used execute
%    faster for this calculation, but it is also more accurate.

end













Thursday, January 21, 2016

Special functions: Incomplete beta functions


All the elements presented are not warranted to be correct or free from defects. 
Please report any errors found to  afstblogs@gmail.com



INCOMPLETE BETA FUNCTIONS

Like the incomplete gamma function, there are four "variations" of the outputs we can get:
  1. Incomplete beta ratio lower tail
  2. Incomplete beta ratio upper tail
  3. Incomplete beta lower tail
  4. Incomplete beta upper tail
In this blog when the nature of the function is not specified we refer to the ratio lower tail that will be the object of the Excel formulas and user defined functions presented


BETA DISTRIBUTION


             

EXCEL FORMULA AND USER DEFINED FUNCTIONS



' Incomplete beta ratio function
Function xlINCBETA(x, a, b)
  With WorksheetFunction
    xlINCBETA = .Beta_Dist(x, a, b, 1)
  End With
End Function 


MATLAB


>> betainc(.5,.5,.2,'lower')
ans =
    0.2682
>> 
betainc(.5,.5,.2,'upper')
ans =
    0.7318






Tuesday, January 19, 2016

Special functions: Incomplete gamma functions


All the elements presented are not warranted to be correct or free from defects. 
Please report any errors found to  afstblogs@gmail.com


The terminology used in literature for the incomplete gamma functions, also called regularized or normalized can be deceiving. In fact there are four "variations" of the function shown in the table below which brings some surprises when comparing results obtained from different applications. For example, when translating code from FORTRAN to BASIC it is quite common to find the former returning upper tails.



In this blog the references to the incomplete gamma function mean the incomplete gamma ratio function.


Excel Gamma distribution pdf




NOTE: It can be found, for example in SPSS manual, the expression of the  Gamma distribution with the parameter β on the numerator. Of course the results obtained are the same being the value of the parameter the inverse between the two expressions. The problem is to use one of the formulas with the β value appropriate for the other.

The incomplete gamma gamma ratio function can be calculated as a particular case of the Gamma cdf distribution (β=1).  We consider only the incomplete gama ratio lower tail since the other can be easily computed from it.  


Excel formulas and user defined functions



' Incomplete gamma ratio function
Function xlGAMMAINC(x, a)
  With WorksheetFunction
    xlGAMMAINC = .Gamma_Dist(x, a, 1, 1)
  End With
End Function



Matlab



>> gammainc(0.5,0.2)
ans =
    0.8788

 >> gammainc(0.5,0.2,'upper')
ans =
    0.1212

>> gammainc(0.5,0.2)+gammainc(0.5,0.2,'upper')
ans =
     1

>> gammainc(0.5,0.2)*gamma(0.2)
ans =
    4.0343




Monday, January 18, 2016

Special functions: Beta


All the elements presented are not warranted to be correct or free from defects. 
Please report any errors found to  afstblogs@gmail.com




Excel formulas and user defined functions




' Beta function
Function xlBETA(a, b)
  With WorksheetFunction
    xlBETA = Exp(.GammaLn_Precise(a) _
      + .GammaLn_Precise(b) _
      - .GammaLn_Precise(a + b))
  End With
End Function

Matlab


>> a=(0.5:5)',b=(0.5:5)',beta(a,b)
a =
                       0.5
                       1.5
                       2.5
                       3.5
                       4.5
b =
                       0.5
                       1.5
                       2.5
                       3.5
                       4.5
ans =
          3.14159265358979
         0.392699081698724
        0.0736310778185108
        0.0153398078788564
       0.00335558297349984


Special functions: Gamma



All the elements presented are not warranted to be correct or free from defects.

Please report any errors found to afstblogs@gmail.com

Use also this address to get the Word and Excel files available concerning this subject

Gamma function
ERRATA
There is a minus sign missing in the formula for negative values as can be found by the values obtained in Maple. The formula presented returns a value with the opposite sign, that is, the value of -GAMMA(-a)
The Gamma(1/2) = sqrt(Pi)
21 DEZ 2017
1. Definition

Although the Gamma function can be defined in the complex plane we will only analyze its computation for real arguments.



Since the function is not defined for a=0, the leftmost a value used in the graph was a=1E+10.  

If the values of the function for any interval ]a,a+1 ] are known, it is possible possible to compute the values for the other intervals using the expression for a<>0.

Γ(a) = (a-1) Γ(a-1) 










2. Extension to negative values of a

However, it is quite common to see the gamma function evaluated for negative values as shown in the next figure that shows that it it is undefined for zero and negative integers.








For a < 0 the function can be evaluated using the following reflection formula:

Excel functions and formulas

The function GAMMALN.PRECISE only works for a > 0.
The formulas if column H use the GAMMALN.PRECISE and the reflection formula to extend the computation to non-positive values. There is no convergence for zero and negative integers. 


DNC - Did not converge

C3: =EXP(GAMMALN.PRECISE(B3))
D3: =EXP(-1)/(GAMMA.DIST(1,B3,1,0))

H3: =IF(F3=0,"DNC",IF(F3<0,IF(F3=INT(F3),"DNC",PI()/
(ABS(F3)*SIN(PI()*ABS(F3))*EXP(GAMMALN.PRECISE(ABS(F3))))), EXP(GAMMALN.PRECISE(F3))))


One alternative formula for column H is to use a very large number instead of DNC in order to facilitate graphing the functions for the a values for which there is no convergence.



However, extending the computation for the entire real field implies non-convergence at the left and the right of zero and negative integers, requiring a solution difficult to implement using formulas in dynamic models.
The graph below was made using six series of paired data sets.










Excel user defined functions





Maple commands and Excel add-in

Maple commands and functions are case sensitive.



L3: =maple("GAMMA(&1)",K3)







'Gamma function
Function xlGAMMA(a)
  With WorksheetFunction
    xlGAMMA = Exp(.GammaLn_Precise(a))
  End With
End Function


Matlab

>> gamma(0.5)
ans =
          1.77245385090552


>> a=(0.5:5),gamma(a)'

>> a=(0.5:5)',GammaF=gamma(a)
a =
                       0.5
                       1.5
                       2.5
                       3.5
                       4.5
GammaF =
          1.77245385090552
         0.886226925452758
          1.32934038817914
          3.32335097044784
          11.6317283965675