|
wolffd@0
|
1 function X2 = chisquared_table(P,v)
|
|
wolffd@0
|
2 %CHISQUARED_TABLE computes the "percentage points" of the
|
|
wolffd@0
|
3 %chi-squared distribution, as in Abramowitz & Stegun Table 26.8
|
|
wolffd@0
|
4 % X2 = CHISQUARED_TABLE( P, v ) returns the value of chi-squared
|
|
wolffd@0
|
5 % corresponding to v degrees of freedom and probability P.
|
|
wolffd@0
|
6 % P is the probability that the sum of squares of v unit-variance
|
|
wolffd@0
|
7 % normally-distributed random variables is <= X2.
|
|
wolffd@0
|
8 % P and v may be matrices of the same size size, or either
|
|
wolffd@0
|
9 % may be a scalar.
|
|
wolffd@0
|
10 %
|
|
wolffd@0
|
11 % e.g., to find the 95% confidence interval for 2 degrees
|
|
wolffd@0
|
12 % of freedom, use CHISQUARED_TABLE( .95, 2 ), yielding 5.99,
|
|
wolffd@0
|
13 % in agreement with Abramowitz & Stegun's Table 26.8
|
|
wolffd@0
|
14 %
|
|
wolffd@0
|
15 % This result can be checked through the function
|
|
wolffd@0
|
16 % CHISQUARED_PROB( 5.99, 2 ), yielding 0.9500
|
|
wolffd@0
|
17 %
|
|
wolffd@0
|
18 % The familiar 1.96-sigma confidence bounds enclosing 95% of
|
|
wolffd@0
|
19 % a 1-D gaussian is found through
|
|
wolffd@0
|
20 % sqrt( CHISQUARED_TABLE( .95, 1 )), yielding 1.96
|
|
wolffd@0
|
21 %
|
|
wolffd@0
|
22 % See also CHISQUARED_PROB
|
|
wolffd@0
|
23 %
|
|
wolffd@0
|
24 %Peter R. Shaw, WHOI
|
|
wolffd@0
|
25 %Leslie Rosenfeld, MBARI
|
|
wolffd@0
|
26
|
|
wolffd@0
|
27 % References: Press et al., Numerical Recipes, Cambridge, 1986;
|
|
wolffd@0
|
28 % Abramowitz & Stegun, Handbook of Mathematical Functions, Dover, 1972.
|
|
wolffd@0
|
29
|
|
wolffd@0
|
30 % Peter R. Shaw, Woods Hole Oceanographic Institution
|
|
wolffd@0
|
31 % Woods Hole, MA 02543 pshaw@whoi.edu
|
|
wolffd@0
|
32 % Leslie Rosenfeld, MBARI
|
|
wolffd@0
|
33 % Last revision: Peter Shaw, Oct 1992: fsolve with version 4
|
|
wolffd@0
|
34
|
|
wolffd@0
|
35 % ** Calls function CHIAUX **
|
|
wolffd@0
|
36 % Computed using the Incomplete Gamma function,
|
|
wolffd@0
|
37 % as given by Press et al. (Recipes) eq. (6.2.17)
|
|
wolffd@0
|
38
|
|
wolffd@0
|
39 [mP,nP]=size(P);
|
|
wolffd@0
|
40 [mv,nv]=size(v);
|
|
wolffd@0
|
41 if mP~=mv | nP~=nv,
|
|
wolffd@0
|
42 if mP==1 & nP==1,
|
|
wolffd@0
|
43 P=P*ones(mv,nv);
|
|
wolffd@0
|
44 elseif mv==1 & nv==1,
|
|
wolffd@0
|
45 v=v*ones(mP,nP);
|
|
wolffd@0
|
46 else
|
|
wolffd@0
|
47 error('P and v must be the same size')
|
|
wolffd@0
|
48 end
|
|
wolffd@0
|
49 end
|
|
wolffd@0
|
50 [m,n]=size(P); X2 = zeros(m,n);
|
|
wolffd@0
|
51 for i=1:m,
|
|
wolffd@0
|
52 for j=1:n,
|
|
wolffd@0
|
53 if v(i,j)<=10,
|
|
wolffd@0
|
54 x0=P(i,j)*v(i,j);
|
|
wolffd@0
|
55 else
|
|
wolffd@0
|
56 x0=v(i,j);
|
|
wolffd@0
|
57 end
|
|
wolffd@0
|
58 % Note: "old" and "new" calls to fsolve may or may not follow
|
|
wolffd@0
|
59 % Matlab version 3.5 -> version 4 (so I'm keeping the old call around...)
|
|
wolffd@0
|
60 % X2(i,j) = fsolve('chiaux',x0,zeros(16,1),[v(i,j),P(i,j)]); %(old call)
|
|
wolffd@0
|
61 X2(i,j) = fsolve('chiaux',x0,zeros(16,1),[],[v(i,j),P(i,j)]);
|
|
wolffd@0
|
62 end
|
|
wolffd@0
|
63 end
|
