Home > wafo > wstats > dist2dsmfun2.m

dist2dsmfun2

PURPOSE ^

Smooths the conditional DIST2D distribution parameters.

SYNOPSIS ^

sphat=dist2dsmfun2(phat,h2,csm,lin,visual,f0)

DESCRIPTION ^

 DIST2DSMFUN2 Smooths the conditional DIST2D distribution parameters.  
   
  CALL:  sphat = dist2dgsmfun2(phat,x2,csm,lin) 
   
     sphat = smoothed parameter structure array evaluated at x2 
     phat  = parameter structure array 
     x2    = evaluation points (default phat.x{1}(:,1))  
     csm   = [csma csmb csmc], smoothing parameter vector which defines 
             the smoothing of parameter A,B and C, respectively  
             0 -> LS straight line 
             1 -> cubic spline interpolant (default [1 1 1]) 
     lin   = [lina linb linc], vector defining the extrapolation of 
             the parameters A,B and C,  respectively (default [1 1 1]) 
             0 : No linear extrapolation outside the range of data 
             1 : Linear extrapolation outside the range of data  
   visual  = 0 : No visual fitting by the user (default) 
             1 : Visual fitting by the user by specifying points using 
                 the mouse and interpolate the points with a spline afterwards. 
  
   See also   dist2dfit smooth

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function sphat=dist2dsmfun2(phat,h2,csm,lin,visual,f0) 
002 %DIST2DSMFUN2 Smooths the conditional DIST2D distribution parameters.  
003 %  
004 % CALL:  sphat = dist2dgsmfun2(phat,x2,csm,lin) 
005 %  
006 %    sphat = smoothed parameter structure array evaluated at x2 
007 %    phat  = parameter structure array 
008 %    x2    = evaluation points (default phat.x{1}(:,1))  
009 %    csm   = [csma csmb csmc], smoothing parameter vector which defines 
010 %            the smoothing of parameter A,B and C, respectively  
011 %            0 -> LS straight line 
012 %            1 -> cubic spline interpolant (default [1 1 1]) 
013 %    lin   = [lina linb linc], vector defining the extrapolation of 
014 %            the parameters A,B and C,  respectively (default [1 1 1]) 
015 %            0 : No linear extrapolation outside the range of data 
016 %            1 : Linear extrapolation outside the range of data  
017 %  visual  = 0 : No visual fitting by the user (default) 
018 %            1 : Visual fitting by the user by specifying points using 
019 %                the mouse and interpolate the points with a spline afterwards. 
020 % 
021 %  See also   dist2dfit smooth  
022  
023 % Tested on: Matlab 5.2 
024 % history 
025 % revised pab July2004   
026 % revised pab 20.01.2004   
027 % revised pab 08.02.2000 
028 %  - added visual 
029 % by  Per A. Brodtkorb 28.10.98 
030  
031  
032 if nargin<2 |isempty(h2), 
033   ind0 = find(~isnan(phat.x{1}(:,1))); 
034   h2=phat.x{1}(ind0,1); 
035 end 
036 h2=h2(:); 
037  
038 if nargin<3|isempty(csm) 
039   csm= []; 
040 end 
041 if (nargin<4) | isempty(lin), 
042     lin=[]; 
043 end 
044 if (nargin<5) | isempty(visual), 
045     visual=0; 
046 end 
047 if (nargin<6)  
048   f0 = []; 
049 end 
050 sphat=phat; 
051  
052 if isfield(phat,'CI') 
053   sphat=rmfield(sphat,'CI'); 
054 end 
055 sphat.date=datestr(now); 
056 pvhsiz=size(phat.x{1}); 
057   
058 n=length(h2(:)); 
059 spvhsiz=[n pvhsiz(2)]; 
060 sphat.x{1}=zeros(spvhsiz); 
061 sphat.x{1}(:,1)=h2; 
062  
063 useSTATS = (1 & isfield(phat,'stats1')); 
064 if useSTATS & ~isempty(phat.stats1) & strncmpi(phat.dist{1},'gamma',2) 
065    CSMA = 1; 
066    CSMB = 1; 
067    linA = 1; 
068    linB = 1; 
069    if (~isempty(csm) & length(csm)>0), CSMA = csm(1); end 
070    if (~isempty(csm) & length(csm)>1), CSMB = csm(2); end 
071    if (~isempty(lin) & length(lin)>0), linA = lin(1); end 
072    if (~isempty(lin) & length(lin)>1), linB = lin(2); end 
073     %smooth on conditional mean and variance 
074     ind1   = find(~isnan(sum([phat.stats1{:}],2))); 
075     hi = phat.stats1{1}(ind1); 
076      
077      
078     mi = phat.stats1{2}(ind1); 
079     vi = phat.stats1{3}(ind1); 
080     Nptsi = phat.stats1{4}(ind1); 
081     dab =  hi.^2./(Nptsi); 
082     %h2 = unique(h2); 
083     Mv = smooth(hi,mi,CSMA,h2,linA,dab); 
084     if any(Mv<=0)  
085       dab1 =  log(mi)+log(mi+dab); 
086        Mv = exp(smooth(hi,log(mi),CSMA,h2,linA,dab1)); 
087     end 
088     Sv = smooth(hi,sqrt(vi),CSMB,h2,linB,dab); 
089     if any(Sv<=0)  
090       dab1 =  log(sqrt(vi))+log(sqrt(vi)+dab); 
091        Sv = exp(smooth(hi,log(sqrt(vi)),CSMB,h2,linB,dab1)); 
092     end 
093     %plot(hi,mi,'.',hi,sqrt(vi),'.',h2,Mv,h2,Sv) 
094     %drawnow 
095     %pause 
096     if 1 
097       Av = Mv.^2./Sv.^2; 
098       Bv = Sv.^2./Mv; 
099     else 
100       Av = smooth(h2,Mv.^2./Sv.^2,0.99,h2); 
101       Bv = smooth(h2,Sv.^2./Mv,0.99,h2); 
102     end 
103     Cv = 0; 
104     sphat.x{1}(:,2) = Av; 
105 else 
106   [sphat.x{1}(:,2), Bv , Cv]=dist2dsmfun(phat,h2,csm,lin); 
107 end 
108   if ~isempty(Bv) 
109     sphat.x{1}(:,3)=Bv;   
110   end 
111  
112   if ~isempty(Cv) & any(Cv~=0), 
113     sphat.x{1}(:,4)=Cv; 
114      
115   end 
116 sphat.csm=csm; 
117 sphat.lin=lin; 
118  
119 if ~visual, return, end 
120  
121 % Visual Fitting of the parameters starts here: 
122 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
123  
124 % this may be as good as any other method when extrapolating the 
125 % parameter outside the observed range. 
126  
127 sphat0 = sphat; 
128  
129 fig0 = gcf; 
130 figure(fig0) 
131 x1 = linspace(0,5); 
132 clf 
133 cmpDistributions(x1,sphat,f0) 
134 figure(fig0+1) 
135 dist2dparamplot(phat,sphat) 
136  
137  
138 for ix = 1:10 
139   figure(fig0+2) 
140   dist2dparamplot(phat,sphat,1) 
141   if (ix>1), 
142     hold on 
143     dist2dparamplot(phat,sphat0,1) 
144     hold off 
145   end 
146   [x,y]=graphicinput; 
147    
148   sphat.visual=~isempty(x); 
149   if ~isempty(x), 
150     sphat.x{1}(:,2)=smooth(x,y,1,sphat.x{1}(:,1),1); 
151     hold on, plot(h2,sphat.x{1}(:,2),'g-'),hold off % plot visual fit 
152     figure(fig0) 
153     cmpDistributions(x1,sphat,f0); 
154     figure(fig0+1) 
155     dist2dparamplot(phat,sphat) 
156   end 
157  
158   if ~isempty(Bv) 
159     figure(fig0+3) 
160     dist2dparamplot(phat,sphat,2) 
161     if (ix>1), 
162       hold on 
163       dist2dparamplot(phat,sphat0,2) 
164       hold off 
165     end 
166     [x,y]=graphicinput; 
167     sphat.visual(2)=~isempty(x); 
168     if ~isempty(x), 
169       sphat.x{1}(:,3)=smooth(x,y,1,sphat.x{1}(:,1),1); 
170       hold on, plot(h2,sphat.x{1}(:,3),'g-'),hold off % plot visual fit 
171       figure(fig0) 
172       cmpDistributions(x1,sphat,f0) 
173       figure(fig0+1) 
174       dist2dparamplot(phat,sphat) 
175     end 
176   end 
177   if ~isempty(Cv) & any(Cv~=0), 
178     if strcmpi(phat.dist{1}(1:2),'ra'), col=2;else,col=3;end 
179     figure(fig0+4) 
180     dist2dparamplot(phat,sphat,col) 
181     if (ix>1), 
182       hold on 
183       dist2dparamplot(phat,sphat0,col) 
184       hold off 
185     end 
186     [x,y]=graphicinput; 
187     sphat.visual(col)=~isempty(x); 
188     if ~isempty(x), 
189       sphat.x{1}(:,col+1)=smooth(x,y,1,sphat.x{1}(:,1),1); 
190       hold on, plot(h2,sphat.x{1}(:,2),'g-'),hold off % plot visual fit 
191       figure(fig0) 
192       cmpDistributions(x1,sphat,f0) 
193       figure(fig0+1) 
194       dist2dparamplot(phat,sphat) 
195     end  
196   end 
197   if ix<10 
198     ButtonName=questdlg('Do you want to visually fit once more?', ... 
199                        'Yes','No'); 
200     if ~strncmpi(ButtonName,'yes',3) 
201       break 
202     end 
203   end 
204 end 
205  
206 return 
207  
208  
209 function [x, y]= graphicinput, 
210 % Local function which reads points selected with the mouse from screen 
211 %  
212   x=[];y=[]; 
213   n=0; 
214   if 0 
215   lim=axis; % find the current scaling of the figure 
216   %answer = input(['Enter scaling for the current plot. (default [' num2str(lim,4) '] )']) 
217   prompt={'Enter scaling for the current plot.'}; 
218   answer=inputdlg(prompt,'Visual fitting',1,{ num2str(lim,4) },'on');     
219   Na=length(answer); 
220   if Na>0, 
221     answer{1}=str2num(answer{1}); 
222     if (isempty(answer{1})|all(answer{1}==lim)),  else, axis(answer{1}), end 
223   end 
224   end 
225   disp('Left mouse button picks points') 
226   disp('Middle mouse button (or z) to zoom in/out') 
227   disp('Right mouse button picks last point') 
228   disp('(if less than 6 points are selected then no fitting is made.)') 
229    
230   leftButton    = 1; 
231   middleButton  = 2; 
232   button = 1; 
233  
234   hold on 
235   while (~isempty(button) & (button==leftButton)), 
236     [xi,yi,button] = ginput(1); 
237     if ~isempty(button) 
238       if (button == leftButton), 
239     plot(xi,yi,'go') 
240     n=n+1; 
241     x(n,1)=xi; 
242     y(n,1)=yi; 
243       elseif( (button == middleButton) | strcmpi(char(button),'z') ) 
244     disp('click the left mouse button to zoom in on the point') 
245     disp('under the mouse.') 
246     disp('Click the right mouse button to zoom out.') 
247     disp('Click and drag to zoom into an area.') 
248     %disp('To resume the identification of points:) 
249     zoom on 
250     button = input('Hit the ENTER key on keyboard when finished zooming.'); 
251     zoom off 
252     button = leftButton; 
253       end 
254     end 
255   end 
256    
257   if n<6, 
258     x=[];y=[]; 
259   end 
260   hold off 
261  
262 return 
263  
264 function cmpDistributions(x1,sphat,f0) 
265  f = dist2dpdf2(x1,x1,sphat); 
266  pdfplot(f); 
267  if ~isempty(f0),  
268    hold on 
269    pdfplot(f0,'r') 
270    hold off 
271  end 
272  return

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