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:
 wtweibcdf Truncated Weibull cumulative distribution function drawnow Flush pending graphics events. mean Average or mean value. num2str Convert number to string. (Fast version) plot Linear plot. polyval Evaluate polynomial.
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