Re: fitting 2d normal distribution

[Francesco Potorti`]

On 23-Oct-2008, Søren Hauberg wrote:
>The attached version extends the
>documentation, and adds a demonstration, plus some minor code changes to
>plot the ellipses if less than two outputs are requested (if one is
>requested a handle is returned).

That's wonderful!  Thanks to all that have contributed to this piece of
software.

>The attached version only plots the ellipse, and its major/minor axes,
>because that's what I prefer. But as I said that's really a matter of
>taste, so there should probably be some way of changing what
>specifically is plotted. 

Yes.  For my application, I think that having the ellipses only is
best.  Also, given the name of the function, its basic usage should be
that of simply drawing an ellipse.  If you agree with this, the default
should be to not drawing anything else.  Also, the help string should
more clearly explain what is the AMAT argument.  Maybe semiaxes could be
a better name than amat?

>I'm just not sure what the best API for doing
>this might be. Suggestions are more than welcome...

I see at least two ways, the first one with two variations:

1) you add an optional fifth argument
1a) You use it as a three-value option meaning
    - ellipse only
    - ellipse with axes
    - ellipse, axes, bbox
1b) You use it as a three-options on-off variables, allowing one to
    choose among the seven possible significant combinations of drawing
    one or more of ellipse, axes, bbox
2) when additional arguments are passed to plot, you check for them and
   look for possible "ellipse", "axes", "bbox" strings followed by "on"
   or "off" and act accordingly

Here is an example of usage when you start from an image showing
distribution of x-y points as shades of gray.  Would it make sense to
make a couple of function from this example, as detailed in the
comments? 

## Generate data as in ellipse's documentation, with 10000 points
X = randn (10000, 2);
scale = 2;
X (:, 1) *= scale;
theta = 0.4;
R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
X = X * R;

## Move them to (x,y)=(100;50)
X = round(10*X);
X(:,1) += 100;
X(:,2) += 50;
x = 0:200; y = 0:100;

## Now create an image showing the frequency of observed values: this is
## the equivalent of a 2d hist function.  Should not we have something
## like that as a function?
X(X<0) = 0;
X(X(:,1)>200,1) = 200;
X(X(:,2)>100,2) = 100;
X += 1;
img = accumarray(fliplr(X),ones(length(X),1),[101 201]);

## Now we have the synthetically generated image that in real life
## should be our starting data and the x and y array describing the
## values.  Let's see how to display it together with the ellipse
## representing its 2nd order statistics.  Should this also be a
## ready-made function?

## Scale the image and compute 1st and 2nd order statistics
img /= sum(img(:));
weightx = sum(img,1);
meanx = weightx*x';
varx = weightx*(x-meanx)'.^2;
weighty = sum(img,2);
meany = y*weighty;
vary = (y-meany).^2*weighty;
[xx yy] = meshgrid(x-meanx,y-meany);
covxy = sum((xx.*yy.*img)(:));

## Now display the image and the superimposed ellipsis
imagesc(x,y,img);  axis("xy", "manual", "equal");
colormap(flipud(gray())); hold on;
ellipse(inv([varx covxy;covxy vary]),2,[meanx meany],100,
        "linewidth",2,"color","black");
hold off

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/