From: redclab on
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
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
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
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
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