001 function [Av , Bv , Cv ]=dist2dsmfun(phat,h2,CSMA,lin_extrap)
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
033
034 PVH = phat.x{1};
035 ind = find(~isnan(sum(PVH(:,:),2)));
036 PVH = PVH(ind,:);
037 CDIST = phat.dist{1};
038 if nargin<2 |isempty(h2),
039 h2=PVH(:,1);
040 end
041
042 CSMB=1;
043 CSMC=1;
044 if nargin<3|isempty(CSMA)
045 CSMA = 1;
046 elseif length(CSMA)==2
047 CSMB=CSMA(2);
048 CSMA=CSMA(1);
049 elseif length(CSMA)>=2
050 CSMB=CSMA(2);
051 CSMC=CSMA(3);
052 CSMA=CSMA(1);
053 end
054 linB=1;
055 linC=1;
056 if (nargin<4) | isempty(lin_extrap),
057 linA=1;
058 elseif length(lin_extrap)==2,
059 linA=lin_extrap(1);
060 linB=lin_extrap(2);
061 elseif length(lin_extrap)>=2,
062 linA=lin_extrap(1);
063 linB=lin_extrap(2);
064 linC=lin_extrap(3);
065 end
066
067
068 useSIG = (1 & isfield(phat,'CI')) ;
069 useSTATS = (useSIG & isfield(phat,'stats1') & ~isempty(phat.stats1));
070 if useSTATS
071 Npts = phat.stats1{end}(ind);
072 end
073 if useSIG
074 useSIG = any(~isnan(phat.CI{1}(:)));
075 end
076 K5=1;
077 if all(PVH(:,1)>=0)
078 da = (PVH(:,1)).^K5.*ones(size(PVH(:,1)));
079 else
080 da = ones(size(PVH(:,1)));
081 end
082 if useSIG,
083 CIa = phat.CI{1}(ind,1:2);
084 if useSTATS
085 NptsA = Npts;
086 end
087 k1 = find(PVH(:,2)<CIa(:,1)| CIa(:,2)< PVH(:,2));
088 if any(k1)
089 CIa(k1,:) = [];
090 PVH(k1,:)=[];
091 ind(k1) = [];
092 da(k1) = [];
093 if useSTATS
094 NptsA(k1) = [];
095 end
096 end
097 tmp = (diff(CIa,1,2)/4).^2;
098
099 if ~isempty(tmp)
100 k0 = find(~isnan(tmp));
101 maxStd = 10;
102 if any(k0)
103 da(k0) = tmp(k0);
104 maxStd = max(maxStd,10*max(tmp(k0)));
105 end
106 k1 = find(isnan(tmp));
107 if any(k1)
108 da(k1) = maxStd;
109 end
110
111
112 da = da+[0; da(1:end-1)]+...
113 [da(2:end);da(end)]+...
114 PVH(:,1).*[0; abs(diff(PVH(:,2)))];;
115 if useSTATS
116 da = da./sqrt(NptsA);
117 end
118 da = PVH(:,1).*da./(PVH(:,2)+eps);
119 end
120 end
121
122 Av=smooth(PVH(:,1),PVH(:,2),CSMA,h2,linA,da);
123 switch lower(CDIST(1:2)),
124 case {'ra','tr','tg','gu','ga','we','tw','gg'},
125 if (any(Av<0)),
126 da2 = abs(log(PVH(:,2)) - log(PVH(:,2)+ da));
127 Av=exp(smooth(PVH(:,1),log(PVH(:,2)),CSMA,h2,linA,da2));
128 end
129 end
130 Bv=[];
131 if ~strcmp(lower(CDIST(1:2)),'ra')
132 if all(PVH(:,1)>=0)
133 db = PVH(:,1).^K5.*ones(size(PVH(:,1)));
134 else
135 db = ones(size(PVH(:,1)));
136 end
137 if useSTATS
138 NptsB = Npts;
139 end
140 if useSIG
141 CIb = phat.CI{1}(ind,3:4);
142 k1 = find(PVH(:,3)<CIb(:,1)| CIb(:,2)< PVH(:,3));
143 if any(k1)
144 CIb(k1,:) = [];
145 PVH(k1,:)=[];
146 ind(k1) = [];
147 db(k1) = [];
148 if useSTATS
149 NptsB(k1) = [];
150 end
151 end
152
153
154 tmp=(diff(CIb,1,2)/4).^2;
155 if ~isempty(tmp)
156 k0 = find(~isnan(tmp));
157 maxStd = 10;
158 if any(k0)
159 db(k0) = tmp(k0);
160 maxStd = max(maxStd,10*max(tmp(k0)));
161 end
162 k1 = find(isnan(tmp));
163 if any(k1)
164 db(k1) = maxStd;
165 end
166 db = (db+[0; db(1:end-1)]+...
167 [db(2:end);db(end)])/3+...
168 PVH(:,1).*[0; abs(diff(PVH(:,3)))];
169 if useSTATS
170 db = db./sqrt(NptsB);
171 end
172 db = PVH(:,1).*db./(PVH(:,3)+eps);
173
174 end
175 end
176
177 Bv=smooth(PVH(:,1),PVH(:,3),CSMB,h2,linB,db);
178 switch CDIST(1:2) ,
179 case {'lo','ga','we','tw','gg'},
180 if (any(Bv<0)),
181 db2 = abs(log(PVH(:,3)) - log(PVH(:,3)+ db));
182 Bv=exp(smooth(PVH(:,1),log(PVH(:,3)),CSMB,h2,linB,db2));
183 end
184 end
185 end
186
187 Cv=0;
188 if (((size(PVH,2)>=4) & (~strcmpi(CDIST(1:2),'ra'))) | ...
189 ((size(PVH,2)>=3) & (strcmpi(CDIST(1:2),'ra')))),
190 Cv=smooth(PVH(:,1),PVH(:,4),CSMC,h2,linC);
191
192 if all(PVH(:,1)>=0)
193 dc = PVH(:,1).^K5.*ones(size(PVH(:,1)));
194 else
195 dc = ones(size(PVH(:,1)));
196 end
197
198 if useSIG
199 CIc = phat.CI{1}(ind,5:6);
200 if useSTATS
201 NptsC = Npts;
202 end
203 k1 = find(PVH(:,4)<CIc(:,1)| CIc(:,2)< PVH(:,4));
204 if any(k1)
205 CIc(k1,:) = [];
206 PVH(k1,:)=[];
207 ind(k1) = [];
208 dc(k1) = [];
209 NptsC(k1) = [];
210 end
211 tmp = (diff(CIc,1,2)/4).^2;
212
213 if ~isempty(tmp)
214 k0 = find(~isnan(tmp));
215 maxStd = 10;
216 if any(k0)
217 dc(k0) = tmp(k0);
218 maxStd = max(maxStd,10*max(tmp(k0)));
219 end
220 k1 = find(isnan(tmp));
221 if any(k1)
222 dc(k1) = maxStd;
223 end
224 dc = dc+[0; dc(1:end-1)]+...
225 [dc(2:end);dc(end)];
226
227 if useSTATS
228 dc = dc./sqrt(NptsC);
229 end
230 dc = PVH(:,1).*dc./(PVH(:,4)+eps);
231 end
232 end
233 switch lower(CDIST(1:2)) ,
234 case {'ra','tg','lo','ga','we'},
235 if (any(Cv<0)),
236 Cv(Cv<0)=0;
237 end
238 case 'gg',
239 if (any(Cv<0)),
240 Cv=exp(smooth(PVH(:,1),log(PVH(:,4)),CSMC,h2,linC,dc));
241 end
242 end
243 end
244 return
245
246
247