cochinv(p,k,v)

function x = cochinv(p,k,v)
%COCHINV Inverse of the Cochran's cumulative distribution function (cdf).
%    COCHINV(P,K,V) returns the inverse of the Cochran's C cumulative 
%    distribution function of K variances and V degrees of freedom at
%    the values in P.
%    [Gives the critical value associated to the Cochran's C statistics 
%    ratio of the largest S^2 to their total (max(S^2)/sum(S^2)), as an
%    alternative procedure to test homocedasticity (Cochran, 1941)].
%
%   Syntax: function x = cochinv(p,k,v) 
%
%     Inputs:
% 	        p - accumulative probability value
% 	        k - number of variances
% 	        v - degrees of freedom for each variance: sample sizes must
%               be equal
%
%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             [e-mail]
%
%  May 27, 2009.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2009). Cochinv: Inverse of 
%  the Cochran's cumulative distribution function. A MATLAB file. [WWW document].
%  URL http://www.mathworks.com/matlabcentral/fileexchange/24283
%
%  References:
% 
%  Cochran, W. G. (1941), The distribution of the largest of a set of
%      estimated variances as a fraction of their total. Annals of
%      Eugenics, 11:47-52.
%  Eisenhart, C., M.W. Hastay and W.A. Wallis. (1947), Selected Techniques
%      of Statistical Analysis for Scientific and Industrial Research and
%      Production and Management Engineering. NY:McGraw-Hill.
%

a = 1-p;

if nargin <  3, 
    error('stats:betainv:TooFewInputs','Requires three input arguments.'); 
end

if numel(p) ~= 1 || p <= 0 || p >= 1
    error('stats:cochinv:BadPvalue','P must be a scalar between 0 and 1.');
end

[errorcode p k v] = distchck(3,p,k,v);

if errorcode > 0
   error('Requires non-scalar arguments to match in size.');
end
 
A = v/2;
B = v*(k-1)/2;
C = a/k;
x = betainv(1-C,A,B);