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.

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. mdist2dcdf Joint 2D CDF due to Plackett mdist2dpdf Joint 2D PDF due to Plackett given as f{x1}*f{x2}*G(x1,x2;Psi). mdist2dstat Mean and variance for the MDIST2D distribution. wlogninv Inverse of the Lognormal distribution function wnorminv Inverse of the Normal distribution function error Display message and abort function. inf Infinity. interp1 1-D interpolation (table lookup) linspace Linearly spaced vector. lower Convert string to lowercase. nan Not-a-Number. num2str Convert number to string. (Fast version)
This function is called by:
 mdist2drnd Random points from a bivariate MDIST2D distribution

## 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 %
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