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
Check if all input arguments are either scalar or of common size. | |
Joint 2D CDF due to Plackett | |
Joint 2D PDF due to Plackett given as f{x1}*f{x2}*G(x1,x2;Psi). | |
Mean and variance for the MDIST2D distribution. | |
Inverse of the Lognormal distribution function | |
Inverse of the Normal distribution function | |
Display message and abort function. | |
Infinity. | |
1-D interpolation (table lookup) | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Not-a-Number. | |
Convert number to string. (Fast version) |
Random points from a bivariate MDIST2D distribution |
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
Comments or corrections to the WAFO group