wolffd@0: function X2 = chisquared_table(P,v)
wolffd@0: %CHISQUARED_TABLE computes the "percentage points" of the
wolffd@0: %chi-squared distribution, as in Abramowitz & Stegun Table 26.8
wolffd@0: %   X2 = CHISQUARED_TABLE( P, v ) returns the value of chi-squared 
wolffd@0: %   corresponding to v degrees of freedom and probability P.
wolffd@0: %   P is the probability that the sum of squares of v unit-variance
wolffd@0: %   normally-distributed random variables is <= X2.
wolffd@0: %   P and v may be matrices of the same size size, or either 
wolffd@0: %   may be a scalar.
wolffd@0: %
wolffd@0: %   e.g., to find the 95% confidence interval for 2 degrees
wolffd@0: %   of freedom, use CHISQUARED_TABLE( .95, 2 ), yielding 5.99,
wolffd@0: %   in agreement with Abramowitz & Stegun's Table 26.8
wolffd@0: %
wolffd@0: %   This result can be checked through the function
wolffd@0: %   CHISQUARED_PROB( 5.99, 2 ), yielding 0.9500
wolffd@0: %
wolffd@0: %   The familiar 1.96-sigma confidence bounds enclosing 95% of
wolffd@0: %   a 1-D gaussian is found through 
wolffd@0: %   sqrt( CHISQUARED_TABLE( .95, 1 )), yielding 1.96
wolffd@0: %
wolffd@0: %   See also CHISQUARED_PROB
wolffd@0: %
wolffd@0: %Peter R. Shaw, WHOI
wolffd@0: %Leslie Rosenfeld, MBARI
wolffd@0: 
wolffd@0: % References: Press et al., Numerical Recipes, Cambridge, 1986;
wolffd@0: % Abramowitz & Stegun, Handbook of Mathematical Functions, Dover, 1972.
wolffd@0: 
wolffd@0: % Peter R. Shaw, Woods Hole Oceanographic Institution
wolffd@0: % Woods Hole, MA 02543  pshaw@whoi.edu
wolffd@0: % Leslie Rosenfeld, MBARI
wolffd@0: % Last revision: Peter Shaw, Oct 1992: fsolve with version 4
wolffd@0: 
wolffd@0: % ** Calls function CHIAUX  **
wolffd@0: % Computed using the Incomplete Gamma function,
wolffd@0: % as given by Press et al. (Recipes) eq. (6.2.17)
wolffd@0: 
wolffd@0: [mP,nP]=size(P);
wolffd@0: [mv,nv]=size(v);
wolffd@0: if mP~=mv | nP~=nv, 
wolffd@0:   if mP==1 & nP==1,
wolffd@0:     P=P*ones(mv,nv);
wolffd@0:   elseif mv==1 & nv==1,
wolffd@0:     v=v*ones(mP,nP);
wolffd@0:   else
wolffd@0:     error('P and v must be the same size')
wolffd@0:   end
wolffd@0: end
wolffd@0: [m,n]=size(P);  X2 = zeros(m,n);
wolffd@0: for i=1:m,
wolffd@0:  for j=1:n,
wolffd@0:   if v(i,j)<=10, 
wolffd@0:    x0=P(i,j)*v(i,j);
wolffd@0:   else
wolffd@0:    x0=v(i,j);
wolffd@0:   end
wolffd@0: % Note: "old" and "new" calls to fsolve may or may not follow 
wolffd@0: % Matlab version 3.5 -> version 4 (so I'm keeping the old call around...)
wolffd@0: %   X2(i,j) = fsolve('chiaux',x0,zeros(16,1),[v(i,j),P(i,j)]); %(old call)
wolffd@0:    X2(i,j) = fsolve('chiaux',x0,zeros(16,1),[],[v(i,j),P(i,j)]);
wolffd@0:  end
wolffd@0: end