Truncated Weibull cumulative distribution function | |
Flush pending graphics events. | |
Average or mean value. | |
Convert number to string. (Fast version) | |
Linear plot. | |
Evaluate polynomial. |
001 function y=wtweibfun(phat,x,F,def,monitor) 002 % WTWEIBFUN Is an internal routine for wtweibfit 003 % 004 005 b =2; c1 = 0; 006 007 a = phat(1); 008 Np = length(phat); 009 if Np>1, b = phat(2);end 010 if Np>2, c1 = phat(3); end 011 012 c = abs(c1); 013 014 N = length(F); 015 %monitor = logical(1); 016 switch def 017 case 1, % fit to sqrt(-log(1-F)) 018 if monitor 019 plot(x,sqrt(-log(1-F)),.... 020 x,sqrt(((x+c)./a).^b-(c/a).^b)); drawnow 021 end 022 023 y=mean((-sqrt(-log(1-F))+... 024 sqrt(((x+c)./a).^b-(c/a).^b)).^2) + N*c*(c1<0); 025 case 2, % fit to (-log(1-F)) 026 if monitor 027 plot(x,(-log(1-F)),... 028 x,(((x+c)./a).^b-(c/a).^b)); drawnow 029 end 030 y=mean((-(-log(1-F))+... 031 (((x+c)./a).^b-(c/a).^b)).^2)+N*c*(c1<0); 032 033 case 3, % fit to (-log(1-F)).^(1/b) 034 if monitor 035 plot(x,(-log(1-F)).^(1/b),x,... 036 (((x+c)./a).^b-(c/a).^b).^(1/b)); drawnow 037 end 038 y=mean((-(-log(1-F)).^(1/b)+... 039 (((x+c)./a).^b-(c/a).^b).^(1/b)).^2)+N*c*(c1<0); 040 case 4, % fit to (-log(1-F)+(c/a).^b).^(1/b) 041 if monitor 042 plot(x,(-log(1-F)+(c/a).^b).^(1/b),x,(x+c)./a); drawnow 043 end 044 y=mean((-(-log(1-F)+(c/a).^b).^(1/b)+(x+c)./a).^2)+N*c*(c1<0); 045 case 5, % fit x/a to ((-log(1-F)+abs(a)).^(1/b)); 046 047 tmp = ((-log(1-F)+abs(a)).^(1/b))-abs(a)^(1/b); 048 p = ([x ]\tmp).'; % Linear LS fit to find 1/a 049 tmp = tmp/p(1); 050 if monitor 051 plot(x,x,x,tmp); drawnow 052 end 053 % Equal weigth on all x: 054 y = (mean(abs((tmp-x)).^(2)))+N*abs(a)*(a<0)+ (b-15)^2*(b>15)/N; 055 case 6, % fit x/a to ((-log(1-F)+abs(a)).^(1/b)); 056 057 cda = abs(a).^(1/b); % = c/a 058 tmp = ((-log(1-F)+abs(a)).^(1/b))-cda; 059 p = ([x ]\tmp).'; % Linear LS fit to find 1/a 060 tmp = tmp/p(1); 061 062 if 0 %monitor 063 plot(x,x,x,tmp); drawnow 064 end 065 066 tmp3 = (-log(1-F)); 067 tmp4 = (((x*p(1)+cda)).^b-abs(a)); 068 069 % fit to (-log(1-F)) 070 % More weight on the tails: The tail is fitted very well 071 y = mean(abs(x-tmp).^(2)+abs(tmp3-tmp4).^(2))+N*abs(a)*(a<0)+(b-6)*(b>10)/N; 072 if monitor 073 plot(x,[x, tmp],x,[tmp3,tmp4]); drawnow 074 end 075 case 7 076 pac=[0.00077598974699 -0.02620368505187 1.28552709525102 -0.73037371897582]; 077 pba=[-0.00641052386506 0.13245900299368 0.45810897559486 -0.38495820627853]; 078 % c = abs((a^1.25-0.4)/1.41)+.2; 079 %c = abs((a^1.25-0.2)/1.45); 080 %a = polyval(pba,b); 081 c = polyval(pac,a); 082 %c = abs((a^1.25-0.57)/1.41); 083 cda = abs(c/a); 084 tmp = (((-log(1-F)+cda^b).^(1/b))-cda)*a; 085 %tmp3 = (-log(1-F)); 086 %tmp4 = ((x+c)/a).^b-cda^b; 087 if monitor 088 plot(x,[x, tmp]); drawnow 089 end 090 y = mean(abs(x-tmp).^(2))+N*abs(a)*(a<=0)+(b-6)*(b>6)/N; 091 case 8 092 tmp = sqrt(-log(1-F)); 093 tmp2 = sqrt(-log(1-wtweibcdf(x,a,b,c))); 094 if monitor 095 plot(x,[ tmp tmp2]); drawnow 096 end 097 y = mean(abs(tmp-tmp2).^(2)); 098 end 099 100 if monitor 101 disp(['err = ' num2str(y,10) ' a b c = ' num2str([a,b,c],4) ]) 102 end 103 104 105 106 107 108 109 110
Comments or corrections to the WAFO group