Last change
on this file since 177 was
37,
checked in by (none), 15 years ago
|
Added original make3d
|
-
Property svn:executable set to
*
|
File size:
1.5 KB
|
Line | |
---|
1 | function y = trigamma(x) |
---|
2 | %TRIGAMMA Trigamma function. |
---|
3 | % TRIGAMMA(X) returns trigamma(x) = d**2 log(gamma(x)) / dx**2 |
---|
4 | % If X is a matrix, returns the trigamma function evaluated at each element. |
---|
5 | |
---|
6 | % Reference: |
---|
7 | % |
---|
8 | % B Schneider, |
---|
9 | % Trigamma Function, |
---|
10 | % Algorithm AS 121, |
---|
11 | % Applied Statistics, |
---|
12 | % Volume 27, Number 1, page 97-99, 1978. |
---|
13 | % |
---|
14 | % From http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f |
---|
15 | |
---|
16 | small = 1e-4; |
---|
17 | large = 8; |
---|
18 | c = pi^2/6; |
---|
19 | c1 = -2.404113806319188570799476; |
---|
20 | b2 = 1/6; |
---|
21 | b4 = -1/30; |
---|
22 | b6 = 1/42; |
---|
23 | b8 = -1/30; |
---|
24 | b10 = 5/66; |
---|
25 | |
---|
26 | % Initialize |
---|
27 | y = zeros(size(x)); |
---|
28 | |
---|
29 | % illegal values |
---|
30 | i = find(isnan(x) | (x == -inf)); |
---|
31 | if ~isempty(i) |
---|
32 | y(i) = nan; |
---|
33 | end |
---|
34 | |
---|
35 | % zero or negative integer |
---|
36 | i = find((x <= 0) & (floor(x)==x)); |
---|
37 | if ~isempty(i) |
---|
38 | y(i) = Inf; |
---|
39 | end |
---|
40 | |
---|
41 | % Negative non-integer |
---|
42 | i = find((x < 0) & (floor(x) ~= x)); |
---|
43 | if ~isempty(i) |
---|
44 | % Use the derivative of the digamma reflection formula: |
---|
45 | % -trigamma(-x) = trigamma(x+1) - (pi*csc(pi*x))^2 |
---|
46 | y(i) = -trigamma(-x(i)+1) + (pi*csc(-pi*x(i))).^2; |
---|
47 | end |
---|
48 | |
---|
49 | % Small value approximation |
---|
50 | i = find(x > 0 & x <= small); |
---|
51 | if ~isempty(i) |
---|
52 | y(i) = 1./(x(i).*x(i)) + c + c1*x(i); |
---|
53 | end |
---|
54 | |
---|
55 | % Reduce to trigamma(x+n) where ( X + N ) >= large. |
---|
56 | while(1) |
---|
57 | i = find(x > small & x < large); |
---|
58 | if isempty(i) |
---|
59 | break |
---|
60 | end |
---|
61 | y(i) = y(i) + 1./(x(i).*x(i)); |
---|
62 | x(i) = x(i) + 1; |
---|
63 | end |
---|
64 | |
---|
65 | % Apply asymptotic formula when X >= large |
---|
66 | i = find(x >= large); |
---|
67 | if ~isempty(i) |
---|
68 | z = 1./(x(i).*x(i)); |
---|
69 | y(i) = y(i) + 0.5*z + (1.0 + z.*(b2 + z.*(b4 + z.*(b6 + z.*(b8 + z.*b10))))) ./ x(i); |
---|
70 | end |
---|
Note: See
TracBrowser
for help on using the repository browser.