| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106 |
- function c = spm_nCr(n,r)
- % Combinatorics: n choose r, nCr
- % FORMAT c = spm_nCr(n,r)
- %
- % n - Number of objects
- % r - Number of objects to choose
- % c - n choose r
- %__________________________________________________________________________
- %
- % spm_nCr returns the number of ways of choosing r objects from a pool
- % of n objects, without replacement, order unimportant. Equivalently:
- % the number of ways of putting r objects into n indistinguishable urns
- % with exclusion. These are the Binomial coefficients of Pascal's
- % triangle:
- % ( n ) n!
- % | | = ---------
- % ( r ) r! (n-r)!
- %
- % n & r must be whole numbers, with n>=r. Non-integer or out-of-range
- % arguments return zero as nCr. Non-scalar n & r must have the same
- % dimensons.
- %
- % Algorithm:
- %--------------------------------------------------------------------------
- % For vary small n, nCr can be computed naively as the ratio of
- % factorials, using gamma(x+1) to return x!. For moderately sized n, n!
- % (& r! &/or (n-r)!) become very large, and naive computation isn't
- % possible. Direct computation is still possible upon noting that the
- % expression cancels down to the product of r fractions:
- %
- % n! n n-1 n-(r-1)
- % --------- = - x --- x ... x -------
- % n! (n-r)! r r-1 1
- %
- % Unfortunately this cunning computation (given at the end of this
- % function) is difficult to vectorise. Therefore we compute the log of
- % nCr as (ln(n!)-ln(r!)-ln((n-r)!), using the log-gamma special
- % function gammaln:
- %
- % nCr = exp( gammaln(n+1) - gammaln(r+1) - gammaln(n-r+1) )
- %
- % The result is rounded to cope with roundoff error for smaller values
- % of n & r. 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) 2008 Wellcome Trust Centre for Neuroimaging
- % Andrew Holmes
- % $Id: spm_nCr.m 5219 2013-01-29 17:07:07Z spm $
- %-Format arguments, note & check sizes
- %--------------------------------------------------------------------------
- ad = [ndims(n);ndims(r)];
- dc = max(ad);
- as = [ [size(n),ones(1,dc-ad(1))];...
- [size(r),ones(1,dc-ad(2))] ];
- sc = max(as);
- xa = prod(as,2)>1;
- if sum(xa)==2 && any(diff(as)), error('non-scalar args must match in size'), end
- %-Computation
- %--------------------------------------------------------------------------
- %-Initialise result to zeros
- c = zeros(sc);
- %-Non zero only where n & r are whole and 0<=r<=n
- Q = find( n==floor(n) & r==floor(r) & n>=r & r>=0 );
- if isempty(Q), return, end
- if xa(1), Qn=Q; else Qn=1; end
- if xa(2), Qr=Q; else Qr=1; end
- %-Compute
- c(Q) = round(exp(gammaln(n(Qn)+1) -gammaln(r(Qr)+1) - gammaln(n(Qn)-r(Qr)+1)));
- %-Return
- %--------------------------------------------------------------------------
- return
- %==========================================================================
- %-Direct computation method: (For interest)
- %==========================================================================
- % The following cunning direct computation is faster than using log
- % gammas, but is rather difficult to vectorise.
- %if r<n/2
- % %-better for small r (less terms), append 1's for 0!
- % c = prod([n:-1:n-r+1,1]./[r:-1:1,1]);
- %else
- % %-better for large r (less terms), append 1's for 0!
- % c = prod([n:-1:r+1,1]./[n-r:-1:1,1]);
- %end
|
