MDIST2DFIT Parameter estimates for MDIST2D data. CALL: Phat = mdist2dfit(x1,x2,dist,alpha,method) phat = structure array containing x = cellarray of distribution parameters x{1} = Phat1 marginal parameters of x1 x{2} = Phat2 marginal parameters of x2 x{3} = Psi the interaction parameter between x1 and x2. CI = Confidence interval of distribution parameters dist = as explained below alpha = ------||---------- method = ------||---------- note = Memorandum string date = date and time of creation x1,x2 = input data dist = list of distributions used in the fitting of X1 and X2, respectively. Options are 'tgumbel', 'gumbel', 'lognormal','rayleigh','weibull','gamma'. alpha = confidence level constant (If not given CI will not be computed) method = fitting method used 'SML' or 'MML' (Simultainous or Marginal Maximum Likelihood method) (default 'SML') Example: phat.x={1 2 2 }; phat.dist={'rayl','rayl'}; [R1,R2] = mdist2drnd(1000,phat); Phat2 = mdist2dfit(R1,R2,{'weibull','rayleigh'},.05,'SML') See also mdist2drnd, mdist2dpdf, mdist2dcdf, mdist2dprb
MDIST log-likelihood function. | |
Parameter estimates for Gamma data. | |
Parameter estimates for Gumbel data. | |
Parameter estimates for Normal data. | |
Inverse of the Normal distribution function | |
Parameter estimates for Rayleigh data. | |
Parameter estimates for Weibull data. | |
Create cell array. | |
Correlation coefficients. | |
String representation of date. | |
Display message and abort function. | |
Multidimensional unconstrained nonlinear minimization (Nelder-Mead). | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Not-a-Number. | |
Current date and time as date number. | |
Convert string matrix to numeric array. | |
Compare strings. | |
MATLAB version number. |
001 function phato=mdist2dfit(V,H,dist,alpha,method) 002 % MDIST2DFIT Parameter estimates for MDIST2D data. 003 % 004 % CALL: Phat = mdist2dfit(x1,x2,dist,alpha,method) 005 % 006 % phat = structure array containing 007 % x = cellarray of distribution parameters 008 % x{1} = Phat1 marginal parameters of x1 009 % x{2} = Phat2 marginal parameters of x2 010 % x{3} = Psi the interaction parameter between x1 and x2. 011 % CI = Confidence interval of distribution parameters 012 % dist = as explained below 013 % alpha = ------||---------- 014 % method = ------||---------- 015 % note = Memorandum string 016 % date = date and time of creation 017 % 018 % x1,x2 = input data 019 % dist = list of distributions used in the fitting of X1 and X2, 020 % respectively. Options are 'tgumbel', 'gumbel', 021 % 'lognormal','rayleigh','weibull','gamma'. 022 % alpha = confidence level constant (If not given CI will not be computed) 023 % method = fitting method used 'SML' or 'MML' (Simultainous or Marginal 024 % Maximum Likelihood method) (default 'SML') 025 % 026 % Example: 027 % phat.x={1 2 2 }; 028 % phat.dist={'rayl','rayl'}; 029 % [R1,R2] = mdist2drnd(1000,phat); 030 % Phat2 = mdist2dfit(R1,R2,{'weibull','rayleigh'},.05,'SML') 031 % 032 % See also mdist2drnd, mdist2dpdf, mdist2dcdf, mdist2dprb 033 034 % References: 035 % Plackett, R. L. (1965) "A class of bivariate distributions." 036 % J. Am. Stat. Assoc. 60. 516-22 037 % [1] Michel K. Ochi, 038 % OCEAN TECHNOLOGY series 6 039 % "OCEAN WAVES, The stochastic approach", Cambridge 040 % 1998 p. 133-134. 041 042 043 % tested on: matlab 5.2 044 % history 045 % Revised pab nov2004 046 % -replaced fmins with fminsearch 047 % revised pab Dec2003 048 % revised pab 8.11.1999 049 % - updated header info 050 % - changed phat from vector to structure 051 % Per A. Brodtkorb 29.01.1999 052 053 054 if (nargin< 3) 055 error('To few input arguments') 056 end 057 CDIST=lower(dist{1}); 058 UDIST=lower(dist{2}); 059 060 if (nargin<5)|isempty(method) 061 method='SML'; % options MML, MOM or SML (simultaneous maximum likelihood) 062 ;% marginal maximum likelihood 063 end 064 if (nargin < 4)|isempty(alpha) 065 alpha = []; 066 else 067 p_int = [alpha/2; 1-alpha/2]; 068 end 069 V = V(:); 070 H = H(:); 071 [n, m] = size(V); 072 [n2, m2] = size(H); 073 if n~=n2 074 error('V and H must have equal size') 075 end 076 077 [ ph, phci]=distfit(H,UDIST,method,alpha); 078 [ pv, pvci]=distfit(V,CDIST,method,alpha); 079 080 rho= corrcoef(V,H); 081 rho=rho(2,1); 082 psi=findPsi(rho); 083 phat=[pv ph psi ]; 084 phatCI=[pvci phci]; 085 switch lower(method(1:3)), 086 case 'sml', %simultaneous MLE 087 pinit=phat; 088 mvrs=version;ix=find(mvrs=='.'); 089 if str2num(mvrs(1:ix(2)-1))>5.2, 090 phat = fminsearch('mdist2dlike',pinit,[],V,H,dist); 091 else 092 phat = fmins('mdist2dlike',pinit,[],[],V,H,dist); 093 end 094 case 'mml',% marginal MLE is already found 095 096 end 097 098 if ~isempty(alpha) % 099 [logL,cov]=mdist2dlike(phat,V,H,dist); 100 sa = diag(cov)'; 101 phatCI = wnorminv(repmat(p_int,1,length(phat)),[phat; phat],[sa;sa]); 102 end 103 phato.x=cell(2,1); 104 phato.CI=cell(1,1); 105 n1=length(pv); 106 [phato.x{1}]=phat(1:n1); 107 [phato.x{2}]=phat(n1+1:end-1); 108 [phato.x{3}]=phat(end); 109 110 [phato.CI{1}]=phatCI; 111 phato.dist=dist; 112 phato.alpha=alpha; 113 phato.method=method; 114 phato.note=[]; 115 phato.date=datestr(now); 116 117 function k=findPsi(rho) 118 % finds psi by a simple iteration scheme 119 rtol=1e-6; %relativ tolerance 120 switch rho 121 case 0, k=1; return 122 case 1, k=inf;return 123 case -1, k=0; return 124 otherwise, 125 if rho<0, 126 kold=0.5; 127 else 128 kold=4; 129 end 130 end 131 kold2=kold-0.001; 132 rho2=getrho(kold2); 133 for ix=1:500, 134 rho1=getrho(kold); 135 k=kold-(rho1-rho).*(kold-kold2)./( rho1-rho2); 136 % disp(['k=' num2str(k) ]), pause 137 if abs(kold-k)<abs(rtol.*k),break 138 elseif k<=0, k=linspace(0,20,1000); 139 [tmp I]=min(abs(getrho(k)-rho)); %linear search 140 k=k(I); 141 break; 142 end 143 kold=k; 144 if ix>=500, disp('could not find k'),end 145 end 146 return 147 148 function r=getrho(psi) 149 r=zeros(size(psi)); 150 k1=find(psi==0); 151 if any(k1), 152 r(k)=-1; 153 end 154 k3=find((psi<1.*psi>0)|psi>1); 155 if any(k3), 156 r(k3)=(psi(k3).^2-1-2.*psi(k3).*log(psi(k3)))./((psi(k3)-1).^2) ; 157 end 158 return 159 160 function [pvhat,pvci]=distfit(tmp,dist,method,alpha) 161 162 if strcmp(lower(method(2:3)),'ml'), 163 switch lower(dist(1:2)), 164 case 'ra', [pvhat tmp2] = wraylfit(tmp,0) ; 165 case 'tg', [pvhat tmp2] = wgumbfit(tmp,0); 166 case 'gu', [pvhat tmp2] = wgumbfit(tmp,0); 167 case 'lo', [pvhat tmp2] = wnormfit(log(tmp),0); 168 case 'ga', [pvhat tmp2] = wgamfit(tmp,0); 169 case 'we', [pvhat tmp2] = wweibfit(tmp,0); 170 otherwise, error('Unknown distribution') 171 end 172 if prod(size(tmp2))~=length(tmp2); 173 sa=diag(tmp2)'; 174 else 175 sa = tmp2(:)'; 176 end 177 p_int=[alpha/2; 1-alpha/2]; 178 pvci = wnorminv(repmat(p_int,1,length(pvhat)),[pvhat; pvhat],[sa;sa]); 179 else % MOM fit 180 pvci=NaN; 181 error('MOM is not implemented for any distribution') 182 switch dist(1:2) 183 case {'tg','gu'} , 184 case 'we', 185 case 'lo', 186 case 'ga', error('MOM is not implemented for Gamma distribution') 187 case 'ra', error('MOM is not implemented for Rayleigh distribution') 188 otherwise , error('Unknown distribution') 189 end 190 end 191 return 192 193
Comments or corrections to the WAFO group