Home > wafo > wstats > mdist2dcinv.m

mdist2dcinv

PURPOSE ^

Inverse of the conditional cdf of X2 given X1.

SYNOPSIS ^

[x, ind]= mdist2dcinv(x1,p,phat);

DESCRIPTION ^

 MDIST2DCINV Inverse of the conditional cdf of X2 given X1.
  
    CALL: x2 =  mdist2dcinv(x1,P,phat)  
 
   x2 = the inverse of mdist2dcdf given x1 and P
   x1 = quantiles
   P  = 0..1
   phat = parameter structure (see mdist2dfit)
 
    The size of X2 is the common size of the input arguments X1 and P. 
    A scalar input functions as a constant matrix of the same size as the
    other inputs.    
 
    MDIST2DCINV uses Newton's method to converge to the solution.
 
  See also   mdist2dcdf,  mdist2dpdf, mdist2dfit, mdist2drnd

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [x, ind]= mdist2dcinv(x1,p,phat);
002 %MDIST2DCINV Inverse of the conditional cdf of X2 given X1.
003 % 
004 %   CALL: x2 =  mdist2dcinv(x1,P,phat)  
005 %
006 %  x2 = the inverse of mdist2dcdf given x1 and P
007 %  x1 = quantiles
008 %  P  = 0..1
009 %  phat = parameter structure (see mdist2dfit)
010 %
011 %   The size of X2 is the common size of the input arguments X1 and P. 
012 %   A scalar input functions as a constant matrix of the same size as the
013 %   other inputs.    
014 %
015 %   MDIST2DCINV uses Newton's method to converge to the solution.
016 %
017 % See also   mdist2dcdf,  mdist2dpdf, mdist2dfit, mdist2drnd
018 
019 %   References:
020 %      Plackett, R. L. (1965) "A class of bivariate distributions."
021 %                                J. Am. Stat. Assoc. 60. 516-22
022 %      [1]  Michel K. Ochi,
023 %       OCEAN TECHNOLOGY series 6
024 %      "OCEAN WAVES, The stochastic approach", Cambridge
025 %      1998 pp. 133-134.
026 
027 
028 %  tested on: matlab 5.2
029 % history
030 % revised pab 8.11.1999
031 %  - updated header info
032 %  - changed phat from vectro to structure
033 % By Per A. Brodtkorb 01.02.99
034 
035 error(nargchk(3,3,nargin))
036 [errorcode  x1 p ] = comnsize(x1,p);
037 
038 if errorcode > 0
039     error('Requires non-scalar arguments to match in size.');
040 end
041 
042 VDIST=lower(phat.dist{1});
043 HDIST= lower(phat.dist{2});
044 
045 %   Initialize X to zero.
046 x = zeros(size(p));
047 %size(x),size(x1)
048 ok=(p > 0  &  p < 1);
049 k = find(~ok);
050 if any(k),
051     tmp  = NaN;
052     x(k) = tmp(ones(size(k))); 
053 end
054 
055 % The inverse cdf of 0 is 0, and the inverse cdf of 1 is inf.  
056 %k0 = find(p == 0);
057 %if any(k0),
058 %    x(k0) = zeros(size(k0)); 
059 %end
060 
061 k1 = find((p == 1| isinf(x1) ));
062 if any(k1), 
063   tmp = Inf;
064     x(k1) = tmp(ones(size(k1))); 
065 end
066 
067 % Newton's Method
068 % Permit no more than count_limit interations.
069 count_limit = 150;
070 count = 0;
071 
072 k = find(~isinf(x1) & ok );
073 
074 pk = p(k);
075 
076 % Supply a starting guess for the iteration.
077 cvar=linspace(eps,max(x1(k)),15); % make sure that the comp. of mn and v do 
078                                 % not consume valuable time
079 [mn v ]=mdist2dstat(phat,1,cvar); %slow
080 mn = interp1(cvar,mn,x1(k),'linear'); %fast
081 v = interp1(cvar,v,x1(k),'linear'); %fast
082 
083 switch 2, %method
084   case 1,%Use a lognormal distribution. 
085   temp = log(v + mn .^ 2); 
086   mu = 2 * log(mn) - 0.5 * temp;
087   sigma = -2 * log(mn) + temp;
088   xk = (wlogninv(pk,mu,sigma.^2));
089   %xk = exp(wnorminv(pk,mu,sigma));
090 case 2, %   Use a normal distribution. probably the fastest choice
091   xk=abs(wnorminv(pk,mn,v));
092   %xk((xk<=0))=1e-3;
093     if any(isnan(xk))
094       disp(['Warning: NaNs   ' num2str(sum(isnan(xk)))])
095     end
096 
097 end
098 %x(k)=xk;
099 %return
100 h = ones(size(pk)); 
101 
102 % Break out of the iteration loop for three reasons:
103 %  1) the last update is very small (compared to x)
104 %  2) the last update is very small (compared to sqrt(eps))
105 %  3) There are more than 100 iterations. This should NEVER happen.
106  
107 k2=find((abs(h) > sqrt(eps)*abs(xk))  &  abs(h) > sqrt(eps));
108 while(any(k2) & count < count_limit), 
109                                  
110     count = count + 1;
111     h(k2)  = (mdist2dcdf(x1(k(k2)),xk(k2),phat,1) - pk(k2)) ./ mdist2dpdf(x1(k(k2)),xk(k2),phat,1);
112    % if any(isnan(h(k2))),      h(isnan(h))=-0.01;    end
113     xnew = xk(k2) - h(k2);
114     % Make sure that xnew>0
115     ix = find(xnew <= 0);
116     if any(ix),
117         xnew(ix) = xk(k2(ix)) / 10;
118         h(k2(ix)) = xk(k2(ix))-xnew(ix);
119     end
120     xk(k2) = xnew;
121      disp(['Iteration  'num2str(count) , '  Number of points left:  ' num2str(length(k2)) ]),
122     %if any(isnan(xk)),  disp(['Warning: values out of range   ' num2str(sum(isnan(xk)))]),   end
123     k2=find((abs(h) > sqrt(eps)*abs(xk))  &  abs(h) > sqrt(eps));   
124 end
125 
126 
127 % Store the converged value in the correct place
128 x(k) = xk;
129 
130 if count == count_limit, 
131     disp('Warning: MDIST2DCINV did not converge.');
132     str = ['The last steps were all less than:  ' num2str(max(abs(h(k2))))];
133     disp(outstr)
134     ind=k(k2);
135   else
136     ind=[];
137 end
138 
139 
140 
141 
142 
143

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003