From: redclab on 23 Mar 2010 06:08 Hi, Im trying to calculate and plot the confidence interval for a custom function I have but I really dont know how I can achieve it, hope you can help me. Here is my code for the regression: dataX=data(:,xcol); dataY=data(:,ycol); k1=1;k2=1;k3=1;k4=1;k5=1;k6=1;k7=1; % Define the exponential function fh=@(x,p) p(1).*(1-(p(2).*exp((-p(3)).*(p(4)./(p(5).^1/3))))).*((p(6).*x).^-p(7)); % define the error function errfh = @(p,x,y) sum((y(:)-fh(x(:),p)).^2); % parameters p0 = [k1 k2 k3 k4 k5 k6 k7]; % search for solution P = fminsearch(errfh,p0,[],dataX,dataY); % % sort the data ord=[squeeze(dataX) fh(dataX,P)]; ord=sortrows(ord); loglog(ord(:,xcol),ord(:,ycol)); With this I get an accurate linear regression (linear because I apply loglog)(I need to have it log scale) but now I need to plot the line with for example 95% confidence and I dont know how to do it. I tried ttest2,normfit..but I get strange values and I dont know how to apply them... Can anyone help me?? thanks in advance
From: Peter Perkins on 23 Mar 2010 09:25 On 3/23/2010 6:08 AM, redclab wrote: > % Define the exponential function > fh=@(x,p) > p(1).*(1-(p(2).*exp((-p(3)).*(p(4)./(p(5).^1/3))))).*((p(6).*x).^-p(7)); Redclab, unless I'm just reading it wrong, this model works out to a constant times x^(-p7). There's no way that you can separately estimate all six of the parameters that go into the leading constant. Reparameterize as p(1)*x^(-p(2)). > With this I get an accurate linear regression (linear because I apply > loglog)(I need to have it log scale) but now I need to plot the line > with for example 95% confidence and I dont know how to do it. I tried > ttest2,normfit..but I get strange values and I dont know how to apply > them... It appears that you have the Statistics Toolbox. I recommend reading this: <http://www.mathworks.com/products/statistics/demos.html?file=/products/demos/shipping/stats/xform2lineardemo.html> and then either using REGRESS on the logged data, or NLINFIT on the original data. Then use NLPREDCI to find confidence intervals.
From: redclab on 24 Mar 2010 09:33 Hi Peter, Thanks for your help... > unless I'm just reading it wrong, this model works out to a constant times x^(-p7). You were not reading wrong, it was me who messed it up with the parameters and added not necessary ones. I read the example you sent on the link (helped me to understand many things) and I rewrote my function. The problem is that after reading the documentation of nlpredci I still dont understand how it works. fh = @(p,x) p(1).*x.^(p(2)); paramEstsLin = [ones(size(x)), log(x)] \ log(y); paramEstsLin(1) = exp(paramEstsLin(1)); xx = linspace(min(x), max(x)); yyLin = fh(paramEstsLin, xx); loglog(x,y,'o', xx,yyLin,'-'); This gives me the regression line Im looking for but now I do the nlinfit to be able to use the nlpredci: [beta,resid,J,Sigma] = nlinfit(x, y, fh, paramEstsLin); [ypred, delta] = nlpredci(fh,x,paramEstsLin,resid,'Covar',Sigma); Delta is supposed to be the confidence interval and i get a vector of 58 elements (the size of my data) that I dont know how to plot. Is that like my 'new y'?? or do I have to add it to the 'y' I already have? or I just shouldnt have got that vector as a ci?? Sorry if my questions are basic but I dont get the concept of it... Thanks, redclab with delta I get a vector of values (of the same size as my data)
From: Peter Perkins on 24 Mar 2010 14:28 On 3/24/2010 9:33 AM, redclab wrote: > [ypred, delta] = nlpredci(fh,x,paramEstsLin,resid,'Covar',Sigma); > > Delta is supposed to be the confidence interval and i get a vector of 58 > elements (the size of my data) that I dont know how to plot. delta is the CI half-width. So [ypred-delta ypred+delta] are the CIs, and the plots usually end up being something like plot(x,y,'bo',x,ypred,'k-',x,ypred-delta,'r:',x,ypred+delta,'r:')
From: redclab on 25 Mar 2010 05:53 Hi Peter, I did what you suggested and it works fine. When I first plotted it, it looked a bit messy. I sorted all the rows of ypred and delta and plotted again and it looks good now. Thank you so much! Greetings, redclab
|
Pages: 1 Prev: Reshape Error Next: "Error closing file" saving structure containing anonymous function |