Home > wafo > wstats > wtweibfun.m

wtweibfun

PURPOSE ^

Is an internal routine for wtweibfit

SYNOPSIS ^

y=wtweibfun(phat,x,F,def,monitor)

DESCRIPTION ^

  WTWEIBFUN Is an internal routine for wtweibfit

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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