001 function sphat=dist2dsmfun2(phat,h2,csm,lin,visual,f0)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
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
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
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
094
095
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
122
123
124
125
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
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
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
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
211
212 x=[];y=[];
213 n=0;
214 if 0
215 lim=axis;
216
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
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