DIST2DFIT Parameter estimates for DIST2D data. CALL: phat = dist2dfit(x1,x2,dist,[res,noverlap C],method,monitor); phat = structure array containing x = cellarray of distribution parameters CI = Confidence interval of distribution parameters dist = as explained below method = ------||---------- note = string date = date and time of creation x1,x2 = input data dist = list of distributions used in the fitting of X1 given X2 and X2, respectively. Options are 'tgumbel', 'gumbel', 'lognormal','rayleigh','trayleigh','weibull','tweibull', 'gamma','ggamma'. res = resolution in the conditional fitting (default range(x2)/12) noverlap = 0,1,2,3..., nr. of groups above and below to overlap in the conditional fitting (default 0) C = defines the location parameter by the following: Cv = min({X1|X2=x}-C,0) if no value is given then Cv=0. method = list containing method used for fitting of X1 given X2 and X2, respectively. Options 'MLE' or 'MOM' (default {'MLE',MLE'}) monitor = 1 if monitor the fitting (default=0) DIST2DFIT fits the following distribution: F{x1|X2=x2}*F{x2}. The parameter(s) of the unconditional distribution of X2, F{x2}, is returned in phat.x{2}. The parameters of the conditional distribution of X1 given X2 are returned in phat.x{1}. The first column in PHAT.X{1} contains the X2 values the parameters in column 2 and 3 are conditioned on. PHAT.CI{1} and PHAT.CI{2} gives 95% confidence interval if MLE are selected for METHOD{1} and METHOD{2}, respectively. Example: R = wraylrnd(2,1000,2); x1 = linspace(0,10)'; phat0=struct('x',[]); phat0.x={[x1,2*ones(size(x1))] 2 }; phat0.dist={'rayl','rayl'}; phat = dist2dfit(R(:,1),R(:,2),{'rayl','rayl'}); sphat = dist2dsmfun2(phat,x1,0); % smooth the parameters dist2dparamplot(phat,sphat); figure(2) % compare fitted model with theoretical f = dist2dpdf2(x1,x1,phat); fs = dist2dpdf2(x1,x1,sphat); fe = dist2dpdf2(x1,x1,phat0); pdfplot(fs); hold on, pdfplot(fe,'k--') plot(R(:,1),R(:,2),'.'),hold off See also dist2drnd, dist2dpdf dist2dcdf dist2dprb
1-dimensional Bin Count | |
Displays a 1D distribution probability plot. | |
Calculates the difference between the maximum and minimum values. | |
Variance | |
Parameter estimates for Gamma data. | |
Parameter estimates for Generalized Gamma data. | |
Parameter estimates for Gumbel data. | |
Plots data on a Gumbel distribution paper. | |
Parameter estimates for Normal data. | |
Inverse of the Normal distribution function | |
Plots data on a Normal distribution paper | |
Parameter estimates for Rayleigh data. | |
Plots data on a Rayleigh distribution paper | |
Parameter estimates for Truncated Rayleigh data. | |
Parameter estimates for truncated Weibull data. | |
Parameter estimates for Weibull data. | |
Plots data on a Weibull distribution paper | |
Control axis scaling and appearance. | |
Create cell array. | |
Create object or return object class. | |
String representation of date. | |
Display message and abort function. | |
Input argument name. | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Average or mean value. | |
Not-a-Number. | |
Current date and time as date number. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Evaluate piecewise polynomial. | |
Print figure or model. Save to disk as image or M-file. | |
Standard deviation. | |
Compare strings ignoring case. | |
Compare first N characters of strings ignoring case. | |
X-axis label. |
setup all global variables of the RECDEMO |
0001 function [phat] =dist2dfit(V,H,dist,res,method,monitor,chat0) 0002 %DIST2DFIT Parameter estimates for DIST2D data. 0003 % 0004 % CALL: phat = dist2dfit(x1,x2,dist,[res,noverlap C],method,monitor); 0005 % 0006 % phat = structure array containing 0007 % x = cellarray of distribution parameters 0008 % CI = Confidence interval of distribution parameters 0009 % dist = as explained below 0010 % method = ------||---------- 0011 % note = string 0012 % date = date and time of creation 0013 % x1,x2 = input data 0014 % dist = list of distributions used in the fitting of X1 given X2 0015 % and X2, respectively. Options are 'tgumbel', 'gumbel', 0016 % 'lognormal','rayleigh','trayleigh','weibull','tweibull', 0017 % 'gamma','ggamma'. 0018 % res = resolution in the conditional fitting (default range(x2)/12) 0019 % noverlap = 0,1,2,3..., nr. of groups above and below to overlap in 0020 % the conditional fitting (default 0) 0021 % C = defines the location parameter by the following: 0022 % Cv = min({X1|X2=x}-C,0) 0023 % if no value is given then Cv=0. 0024 % method = list containing method used for fitting of 0025 % X1 given X2 and X2, respectively. 0026 % Options 'MLE' or 'MOM' (default {'MLE',MLE'}) 0027 % monitor = 1 if monitor the fitting (default=0) 0028 % 0029 % DIST2DFIT fits the following distribution: F{x1|X2=x2}*F{x2}. 0030 % The parameter(s) of the unconditional distribution of X2, 0031 % F{x2}, is returned in phat.x{2}. The parameters of the conditional 0032 % distribution of X1 given X2 are returned in phat.x{1}. The first column 0033 % in PHAT.X{1} contains the X2 values the parameters in column 2 and 3 are 0034 % conditioned on. PHAT.CI{1} and PHAT.CI{2} gives 95% confidence 0035 % interval if MLE are selected for METHOD{1} and METHOD{2}, respectively. 0036 % 0037 % Example: 0038 % R = wraylrnd(2,1000,2); x1 = linspace(0,10)'; 0039 % phat0=struct('x',[]); phat0.x={[x1,2*ones(size(x1))] 2 }; 0040 % phat0.dist={'rayl','rayl'}; 0041 % phat = dist2dfit(R(:,1),R(:,2),{'rayl','rayl'}); 0042 % sphat = dist2dsmfun2(phat,x1,0); % smooth the parameters 0043 % dist2dparamplot(phat,sphat); 0044 % figure(2) % compare fitted model with theoretical 0045 % f = dist2dpdf2(x1,x1,phat); 0046 % fs = dist2dpdf2(x1,x1,sphat); 0047 % fe = dist2dpdf2(x1,x1,phat0); 0048 % pdfplot(fs); hold on, 0049 % pdfplot(fe,'k--') 0050 % plot(R(:,1),R(:,2),'.'),hold off 0051 % 0052 % See also dist2drnd, dist2dpdf dist2dcdf dist2dprb 0053 0054 0055 % tested on: matlab 5.2 0056 % history: 0057 % revised pab Sept 2005 0058 % -fixed abug: [a{1}] = 1 is not allowed anymore!! replaced with a{1} = 1; 0059 % revised pab Jul2004 0060 % Added secret option of C0 to wggamfit 0061 % revised pab 03.12.2000 0062 % -added truncated weibull and rayleigh 0063 % revised pab 12.11.2000 0064 % - added wggampdf option 0065 % - For monitor> 0: removed normplot ...etc and replaced it with a call 0066 % to empdistr instead 0067 % Per A. Brodtkorb 20.10.1998 0068 0069 error(nargchk(3,7,nargin)) 0070 ptime = 1; %pause length if monitor=1 0071 printflag = 0;%print if monitor=1 0072 0073 [HI I] = sort(H(:));%sorting H 0074 VI = V(I); 0075 0076 if (nargin< 3)|isempty(dist), 0077 dist={'tgumbel','rayleigh'}; 0078 else 0079 nd=length(dist); 0080 if nd==1, 0081 dist={dist{1},'rayleigh'}; 0082 end 0083 end 0084 0085 if (nargin<5)|isempty(method) 0086 method={'MLE','MLE'}; % options MLE or MOM 0087 else 0088 nd=length(method); 0089 if nd==1, 0090 method={method{1}, 'MLE'}; 0091 end 0092 end 0093 0094 if (nargin<6)|isempty(monitor) 0095 monitor=0; %if 1 monitor the fitting 0096 end 0097 %monitor=1; 0098 0099 0100 noverlap=0;C=[]; 0101 if nargin<4|isempty(res) 0102 res=range(H(:))/12; % resolution 0103 else 0104 if length(res)>1 0105 noverlap=res(2); 0106 end 0107 if length(res)>2 0108 C=res(3); 0109 end 0110 res=res(1); 0111 end 0112 0113 %N=length(HI); 0114 0115 grp=floor(HI/res)+1; % dividing the data into groups 0116 0117 if monitor 0118 [len,bin] = bincount(grp); 0119 utp = 1/(2*max(len)); 0120 ax = [0 max(VI),utp,1-utp]; 0121 name1=inputname(1); 0122 name2=inputname(2); 0123 if isempty(name1), name1 = 'x1'; end 0124 if isempty(name2), name2 = 'x2'; end 0125 0126 end 0127 [ phhat, phci]=distfit(HI,dist{2},method{2},monitor);%unconditional fitting 0128 if monitor 0129 xlabel(name2) 0130 pause(ptime) 0131 if printflag, print -Pmhlaser ; end 0132 end 0133 0134 Ngrp=grp(end); % # groups 0135 0136 if strcmpi(dist{1}(1:2),'ra'), 0137 npar=1; 0138 elseif strcmpi(dist{1}(1:2),'gg')|strcmpi(dist{1}(1:2),'tw'), 0139 npar = 3; 0140 else 0141 npar=2; 0142 end 0143 %C = input('Specify location parameter (default none)'); 0144 if isempty(C), 0145 pvhat=zeros(Ngrp,1+npar); 0146 else 0147 Cvold=0; 0148 pvhat=zeros(Ngrp,2+npar); 0149 end 0150 0151 pvci=zeros(Ngrp,2*npar); 0152 0153 pvhat(:,1)=linspace(res/2, (Ngrp-0.5)*res, Ngrp)'; 0154 if nargin<7 0155 chat0 = []; 0156 end 0157 0158 Nmin = min(max(6,Ngrp),25); % Minimum number of data in groups 0159 m = zeros(Ngrp,1); 0160 v = m; 0161 Ni = m; 0162 for ix=1:Ngrp, 0163 ind = (grp==ix); 0164 tmp = VI(ind); 0165 Ni(ix) = length(tmp); 0166 if length(tmp)>max(4,0),% if less than 4 observations in the group 0167 m(ix)=mean(tmp); % mean of data in group ix 0168 v(ix)=std(tmp).^2; % variance of data in group ix 0169 else 0170 m(ix)=NaN; 0171 v(ix)=NaN; 0172 end 0173 %if 1 | ((ix-1)*res>2) 0174 % chat0 = 1.5; 0175 %end 0176 switch class(chat0) 0177 case 'inline' 0178 chat00 = chat0(pvhat(ix,1)); 0179 case 'struct' 0180 chat00 = ppval(chat0,pvhat(ix,1)); 0181 otherwise 0182 chat00 = chat0; 0183 end 0184 0185 if ix>=Ngrp-noverlap+1 | ix<=noverlap 0186 for iy=1:min(min(noverlap,ix-1),Ngrp-ix) 0187 kup = (grp==ix+iy); 0188 if any(kup) 0189 ind = (ind | (grp==ix-iy) | kup); 0190 end 0191 end 0192 end 0193 0194 tmp=VI(ind);%find data in group number ix 0195 0196 if length(tmp)<Nmin,% if less than Nmin observations in the group 0197 %grp(grp==ix)=ix-1; %merge the data with group below 0198 %if (ix>3)&~isempty(tmp), grp(grp==ix-2)=ix-1;end % also merge the 0199 % data in grp ix-2 with grp ix-1 0200 pvhat(ix,:)=NaN; 0201 pvci(ix,:)=NaN; 0202 else 0203 0204 if ~isempty(C) 0205 %disp('Cv') 0206 Cv=max(min(tmp)-C,0); 0207 0208 Cv=min(Cv,Cvold+C/30); 0209 tmp=tmp-Cv; 0210 Cvold=Cv; 0211 pvhat(ix,end)=Cv; 0212 end 0213 % conditional fitting 0214 [ pvhat(ix, 2:npar+1),... 0215 pvci(ix,1:2*npar)]=distfit(tmp(:),dist{1},method{1},monitor,chat00); 0216 if monitor 0217 axis(ax) 0218 xlabel([name1 ' | ' name2 '=' num2str(pvhat(ix,1)) ]) 0219 pause(ptime) 0220 if printflag, print -Pmhlaser ; end %print -Pmhlaser 0221 end 0222 end 0223 end 0224 0225 %phat=struct('x',[],'CI',[],'dist',dist,'method',method); 0226 phat.x=cell(2,1); 0227 phat.CI=cell(2,1); 0228 phat.x{1}=pvhat; 0229 phat.x{2}=phhat; 0230 phat.CI{1}=pvci; 0231 phat.CI{2}=phci; 0232 phat.dist=dist; 0233 phat.method=method; 0234 phat.note=[]; 0235 phat.date=datestr(now); 0236 phat.res=res; 0237 phat.noverlap=noverlap; 0238 phat.C=C; 0239 phat.visual=[]; 0240 phat.csm=[]; 0241 phat.lin=[]; 0242 phat.stats1{1} = pvhat(:,1); % x2 0243 phat.stats1{2} = m; %conditional mean given x2 0244 phat.stats1{3} = v; % conditional variance given x2 0245 phat.stats1{4} = Ni;% conditional number of samples given x2 0246 %pvhat(end,:)=[]; 0247 %pvhat(1,:)=[]; 0248 %pvci(end,:)=[]; 0249 %pvci(1,:)=[]; 0250 function [pvhat,pvci]=distfit(tmp,dist,method,monitor,chat) 0251 if nargin>4 & ~isempty(chat) & isfinite(chat) 0252 chat0 = chat; 0253 else 0254 chat0 = []; 0255 end 0256 if strncmpi(method,'mle',3), 0257 switch lower(dist(1:2)), 0258 case 'tr', [pvhat, tmp2] = wtraylfit(tmp,monitor) ; 0259 case 'ra', [pvhat, tmp2] = wraylfit(tmp,monitor) ; 0260 case 'tg', [pvhat, tmp2] = wgumbfit(tmp,monitor); 0261 case 'gu', [pvhat, tmp2] = wgumbfit(tmp,monitor); 0262 case 'lo', [pvhat, tmp2] = wnormfit(log(tmp),monitor); 0263 case 'ga', [pvhat, tmp2] = wgamfit(tmp,monitor); 0264 case 'gg', [pvhat, tmp2] = wggamfit(tmp,monitor,chat0); 0265 %if any(pvhat<0),pvhat(1:3) =NaN; end 0266 case 'we', [pvhat, tmp2] = wweibfit(tmp,monitor); 0267 case 'tw', [pvhat, tmp2] = wtweibfit(tmp,monitor); 0268 %[tmp tmp2]=gumbfit(log(tmp)); 0269 %pvhat=[exp(tmp(2)) 1/tmp(1) ] 0270 otherwise, error('Unknown distribution') 0271 end 0272 if any(isnan(pvhat)) 0273 pvci=[pvhat pvhat]; 0274 else 0275 tz=size(tmp2); 0276 alpha=0.05; % 95% confidense interval 0277 pint = [alpha/2; 1-alpha/2]; 0278 if tz(1)==tz(2), 0279 sa = diag(tmp2)'; 0280 else 0281 sa = tmp2(:)'; 0282 end 0283 nt=length(sa); 0284 %pvhat 0285 pvci = wnorminv(repmat(pint,[1 nt]),[pvhat;pvhat],[sa;sa]); 0286 pvci=pvci(:).'; 0287 end 0288 if 0 %monitor 0289 switch lower(dist(1:2)), 0290 case 'ra',wraylplot(tmp); 0291 case 'gu',wgumbplot(tmp); 0292 case 'tg',wgumbplot(tmp); 0293 case 'lo',wnormplot(log(tmp)); 0294 case 'ga',distplot(tmp,'wgampdf'); 0295 case 'gg',distplot(tmp,'wggampdf'); 0296 case 'we',wweibplot(tmp); 0297 otherwise, error('Unknown distribution') 0298 end 0299 end 0300 0301 else % MOM fit 0302 pvci=NaN; 0303 switch lower(dist(1:2)) 0304 case {'tg','gu'} , pvhat =wgumbplot(tmp); 0305 0306 case 'we', pvhat = wweibplot(tmp); 0307 0308 case 'lo', pvhat = wnormplot(log(tmp)); 0309 0310 case 'ga', ma=mean(tmp);sa=var(tmp); 0311 pvhat=[ma^2/sa sa/ma]; 0312 0313 case 'ra', error('MOM is not implemented for Rayleigh distribution') 0314 otherwise , error('Unknown distribution') 0315 end 0316 pvci = repmat(nan,1,2*length(pvhat)); 0317 0318 if monitor 0319 switch lower(dist(1:2)), 0320 case {'gu' 'tg'}, wgumbplot(tmp); 0321 case 'lo', wnormplot(log(tmp)); 0322 case 'ga', distplot(tmp,'wgampdf'); 0323 case 'we', wweibplot(tmp); 0324 otherwise, error('Unknown distribution') 0325 end 0326 end 0327 end 0328 return 0329 0330 0331 0332 0333 0334
Comments or corrections to the WAFO group