-
Notifications
You must be signed in to change notification settings - Fork 8
/
spm_Ipdf.m
130 lines (119 loc) · 4.63 KB
/
spm_Ipdf.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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
function f = spm_Ipdf(x,n,p)
% Probability Distribution Function of Binomial distribution
% FORMAT f = spm_Ipdf(x,n,p)
%
% x - ordinates
% n - Binomial n
% p - Binomial p [Defaults to 0.5]
% f - PDF
%__________________________________________________________________________
%
% spm_Ipdf returns the Probability (Distribution) Function (PDF) for
% the Binomial family of distributions.
%
% Definition:
%--------------------------------------------------------------------------
% The Bin(n,p) distribution is the distribution of the number of
% successes from n identical independent Bernoulli trials each with
% success probability p. If random variable X is the number of
% successes from such a set of Bernoulli trials, then the PDF f(x) is
% Pr{X=x}, defined for whole n (i.e. non-negative integer n) and
% p in [0,1], given by: (See Evans et al., Ch6)
%
% { nCx * p^x * (1-p)^(n-x) for x=0,1,...,n
% f(x) = |
% { 0 otherwise
%
% where nCx is the Binomial coefficient "n-choose-x", given by n!/(x!(n-x)!).
%
% Algorithm:
%--------------------------------------------------------------------------
% For vary small n, nCx can be computed naively as the ratio of
% factorials, using gamma(n+1) to return n!. For moderately sized n, n!
% (& x! &/or (n-x)!) become very large, and naive computation isn't
% possible. Direct computation is always possible upon noting that the
% expression cancels down to the product of x fractions:
%
% n! n n-1 n-(x-1)
% --------- = - x --- x ... x -------
% n! (n-x)! x x-1 1
%
% Unfortunately this cunning computation (given at the end of this
% function) is difficult to vectorise. Therefore we compute via the log
% of nCx as e^(ln(n!)-ln(x!)-ln((n-x)!), using the special function
% gammaln:
%
% nCx = exp( gammaln(n+1) - gammaln(x+1) - gammaln(n-x+1) )
%
% The result is rounded to cope with roundoff error for smaller values
% of n & x. See Press et al., Sec6.1 for further details.
%
% 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
%
%__________________________________________________________________________
% Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging
% Andrew Holmes
% $Id: spm_Ipdf.m 4182 2011-02-01 12:29:09Z guillaume $
%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<3, p=0.5; end
if nargin<2, error('Insufficient arguments'), end
ad = [ndims(x);ndims(n);ndims(p)];
rd = max(ad);
as = [[size(x),ones(1,rd-ad(1))];...
[size(n),ones(1,rd-ad(2))];...
[size(p),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
error('non-scalar args must match in size');
end
%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
f = zeros(rs);
%-Only defined for whole n, and for p in [0,1]. Return NaN if undefined.
md = ( ones(size(x)) & n==floor(n) & n>=0 & p>=0 & p<=1 );
if any(~md(:))
f(~md) = NaN;
warning('Returning NaN for out of range arguments');
end
%-Non-zero only where defined and x is whole with 0<=x<=n
Q = find( md & x==floor(x) & n>=x & x>=0 );
if isempty(Q), return, end
if xa(1), Qx=Q; else Qx=1; end
if xa(2), Qn=Q; else Qn=1; end
if xa(3), Qp=Q; else Qp=1; end
%-Compute
f(Q) = round(exp(gammaln(n(Qn)+1) -gammaln(x(Qx)+1) - gammaln(n(Qn)-x(Qx)+1)));
f(Q) = f(Q).* p(Qp).^x(Qx) .* (1-p(Qp)).^(n(Qn)-x(Qx));
%-Return
%--------------------------------------------------------------------------
return
%==========================================================================
%-Direct computation method: (For interest)
%==========================================================================
% The following cunning direct computation is faster than using log
% gammas, but is rather difficult to vectorise.
%q=1-p;
%if r<n/2
% %-better for small r (less terms) / small p (smaller numbers)
% %----------------------------------------------------------------------
% f=prod([[n:-1:n-r+1]*p,1]./[r:-1:1,1])*q^(n-r);
%else
% %-better for large r (less terms) / small q (smaller numbers)
% %----------------------------------------------------------------------
% f=prod([[n:-1:r+1]*q,1]./[n-r:-1:1,1])*p^r;
%end