mdist2dcinv

 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

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