-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathspm_Tpdf.m
80 lines (73 loc) · 2.59 KB
/
spm_Tpdf.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function f = spm_Tpdf(x,v)
% Probability Density Function (PDF) of Students t distribution
% FORMAT f = spm_Tpdf(x,v)
%
% x - t-ordinates
% v - degrees of freedom (v>0, non-integer d.f. accepted)
% f - PDF of t-distribution with v degrees of freedom (df) at points t
%_______________________________________________________________________
%
% spm_Tpdf implements the Probability Density Function of Students
% t-distributions.
%
% Definition:
%-----------------------------------------------------------------------
% The Student's t-distribution with v degrees of freedom is defined for
% positive integer v and x in (-Inf,Inf), and has Probability Distribution
% Function (PDF) f(x) given by: (See Evans et al., Ch37)
%
% gamma((v+1)/2)
% f(x) = ----------------------- * (1 + x^2/v) ^ -((v+1)/2
% sqrt(pi*v) * gamma(v/2)
%
% This implementation is not restricted to whole (positive integer) df
% v, rather it will compute for any df v>0.
%
% Algorithm:
%-----------------------------------------------------------------------
% Direct computation using the beta function for
% sqrt(pi)*gamma(v/2) / gamma((v+1)/2) = beta(v/2,1/2)
%
% References:
%-----------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
% "Statistical Distributions"
% 2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
% "Handbook of Mathematical Functions"
% US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
% "Numerical Recipes in C"
% Cambridge
%
%_______________________________________________________________________
% @(#)spm_Tpdf.m 2.2 Andrew Holmes 99/04/26
%-Format arguments, note & check sizes
%-----------------------------------------------------------------------
if nargin<2, error('Insufficient arguments'), end
ad = [ndims(x);ndims(v)];
rd = max(ad);
as = [ [size(x),ones(1,rd-ad(1))];...
[size(v),ones(1,rd-ad(2))]; ];
rs = max(as);
xa = prod(as,2)>1;
if all(xa) & any(diff(as(xa,:)))
error('non-scalar args must match in size'), end
%-Computation
%-----------------------------------------------------------------------
%-Initialise result to zeros
f = zeros(rs);
%-Only defined for v>0. Return NaN if undefined.
md = ( ones(size(x)) & v>0 );
if any(~md(:)), f(~md) = NaN;
warning('Returning NaN for out of range arguments'), end
%-Compute where defined
Q = find( md );
if isempty(Q), return, end
if xa(1), Qx=Q; else Qx=1; end
if xa(2), Qv=Q; else Qv=1; end
%-Compute
f(Q) = ( (1+x(Qx).^2./v(Qv)).^(-(v(Qv)+1)/2) ) ./ ...
(sqrt(v(Qv)).*beta(v(Qv)/2,1/2));