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