[GRASS-dev] Compute mahalanobis distance using Scipy

[Glynn Clements]

Paulo van Breugel wrote:

> > The second version:
> >
> > > >     VI = np.linalg.inv(covar)
> > > >     delta = dat_ref - stat_mean[:,None,None]
> > > >     m = np.sum(np.sum(delta[:,:,:,None] * delta[:,:,None,:] * VI[None,None,:,:],axis=-1),axis=-1)
> > > >     stat_mah = np.sqrt(m)
> >
> > should work with delta being a 3-D array.
> >
> 
> Yes, but when running the lines above, I am getting the same error message:
> 
> Traceback (most recent call last):
>   File "<stdin>", line 1, in <module>
> ValueError: operands could not be broadcast together with shapes
> (3,77,78,78) (1,1,3,3)

Well, it's not quite the same; both operands now have the same number
of dimensions, but they're not compatible (essentially, for each
dimension either both values should be the same or one of them should
be one).

> >From what I can see from broadcasting rules, there is mismatch in the last
> and fore-last dimensions of VI[None,None,:,:] compared to the other two?
> 
> delta[:,:,:,None].shape
> (3, 77, 78, 1)
> 
> delta[:,:,None,:].shape
> (3, 77, 1, 78)
> 
> VI[None,None,:,:].shape
> (1, 1, 3, 3)

The problem was that delta is 3x77x78 whereas the code assumed that it
would be 77x78x3.

> I am really have to get a better grasp of working with these arrays, but in
> any case, after a bit trial and error, I came up with the following to
> compute m.
> 
>  m = np.sum(np.sum(delta[None,:,:,:] * VI[:,:,None,None],axis=1) * delta,
> axis=0)
> 
> 
> This does, in my limited testing, gives the same result as using the loop
> with
> 
> pow(distance.mahalanobis(cell_ref, stat_mean, VI),2).
> 
> 
> to be sure, does the above makes sense to you?

I believe that's correct.

-- 
Glynn Clements <glynn at gclements.plus.com>