Matlab regressions

From LiteratePrograms
Jump to: navigation, search

This program is under development.
Please help to debug it. When debugging
is complete, remove the {{develop}} tag.

A lot of users use Matlab to perform regressions. Let's give some examples.


[edit] Linear regression

You can use the built-in function regress

[B, Bint, R, Rint, Stats] = regress(y, x)

[edit] Apparently nonlinear regressions

[edit] Simple example

[edit] Setting the question

You can use an question asked on the Matlab newsgroup (Help to curve fitting):

I have some data (x,y)and one equation like y=(a*x)/(1+ax) and I need to plot the data and the curve [...]
When you have such a parametrized equation:

(1)\;\; y = \frac{ax}{1+ax}

you can rewrite it in a

(2)\;\; a = \frac{y}{x-x\cdot y}

so if you have some (x,y) distribution, you can plot them and use \mathbf{E}(a) (with a define as in (2)).

result of the trivial estimation script

First generate a test set:

  • f is defined as an anonymous function
  • the real value of a is arbitrary chosen
  • a noise coefficient s is chosen
  • the domain of x is [0,10], the chosen resolution for x is 0.05
  • y is the theoretical value of y
  • ny the noisy (observed) values for y
<<generate test set>>=
f  = @(a,x)(a*x./(1+a*x))
a  = sqrt(2)
s  = .05;
x  = (0:.1:10)';

y  = f(a,x);
ny = y + randn(size(y))*s;

figure('Color',[0.9412 0.9412 0.9412 ]);
h1 = subplot(1,2,1);

[edit] Simple answer

Now we can get the implied values for a:

<<a implied values>>=
g  = @(x,y)(y./(x-x.*y))

ai = g(x,ny);

h2 = subplot(1,2,2);

And a trivial estimator for a is the empirical expectation (the mean) of the implied a:

<<trivial estimation>>=
ah = mean(ai(~isinf(ai)));

title(sprintf('$\\hat a$ = %5.5f ; real a = %5.5f', ah, a), 'interpreter', 'latex');

hold on
plot(x, f(ah, x), ':r', 'linewidth', 2);
hold off
legend({'theoretical', 'observed', 'estimation'});

Let's produce the whole code:


generate test set

a implied values

trivial estimation

(see the upper figure for the result).

But it's a not so simple answer, because from a statistical point of view, you have to master some of the properties of the estimator:

  • is it efficient?
  • is it biaised?
  • etc

[edit] External links

Download code