А, ну тогда все ясно. 
Кроме формулы в статистике... 
Если кому интересно, то вот реализация:
мнения и найденные ошибки приветствуются 
| Код |
{ Здесь некоторая избыточность в переменных
это только для этапа разработки
}
type TMatrix = array of array of double;
TVector = array of double;
function Mahalanobis2(const et, eksp : TMatrix):double; // эталон, эксперимент
var
etl, exp, cov1, cov2, s : TMatrix; // эталон и эксперимент
x1, x2, tmpx, tmpy : TVector; // векторы средних
ai,aj, bi,bj, ci,cj, i,j : integer;
colcount, rowcountexp, rowcountetl : integer;
begin
ErrorMatrix:=false;
result:=-1;//error
rowcountetl:=RowMatrix(et);// количество строк матрицы
rowcountexp := RowMatrix(eksp);// bi
aj:=ColMatrix(et); bj:=ColMatrix(eksp); // количество столбцов матрицы
if aj<>bj then
begin
Messages('Mahalanobis: Количество столбцов у матриц не эквивалентно!');
result:=-1;
ErrorMatrix:=True;
Exit;
end;
etl:=CloneMatrix(et); // создание точной копии матрицы
exp:=CloneMatrix(eksp);
if ErrorMatrix then Exit; //произошла ошибка при клонировании
colcount:=ColMatrix(exp); // просто для удобства
x1:=CreateVector(colcount);
x2:=CreateVector(colcount);
tmpx:=CreateVector(rowcountetl);
tmpy:=CreateVector(rowcountexp);
for i:=0 to colcount-1 do
begin
for j:=0 to rowcountetl-1 do //вычисляем вектора средних
tmpx[j]:=etl[j,i];
for j:=0 to rowcountexp-1 do
tmpy[j]:=exp[j,i];
x1[i]:=mean(tmpx);
x2[i]:=mean(tmpy);
for j:=0 to rowcountetl-1 do //центрируем матрицы
etl[j,i]:=etl[j,i]-x1[i];
for j:=0 to rowcountexp-1 do
exp[j,i]:=exp[j,i]-x2[i];
end;
SetLength(tmpx,0);
SetLength(tmpy,0);
cov1:=MultMatrix(TransMatrix(etl),etl); // умножение и транспонирование матриц
cov2:=MultMatrix(TransMatrix(exp),exp);
s:=CMultMatrix((1/(rowcountetl+rowcountexp-2)),SumMatrix(cov1,cov2)); // умножение на скаляр и суммирование матриц
SetLength(cov1,0);
SetLength(cov2,0);
cov1:=ReversMatrix(s);// обращаем матрицу
Setlength(s,0);
s:=cov1; // просто для удобства
tmpx:=SubVector(x1,x2);// разница векторов
ci:=ColMatrix(s);
cj:=RowMatrix(s);
SetLength(tmpy, cj); // первое умножение матрицы на вектор
for i:=0 to ci-1 do
begin
tmpy[i]:=0;
for j:=0 to cj-1 do
tmpy[i]:=tmpy[i]+s[i,j]*tmpx[j];
end;
result:=MultSVector(tmpy,tmpx); // умножение вектора на вектор
SetLength(s,0);
SetLength(tmpx,0);
setLength(tmpy,0);
setLength(x1,0);
SetLength(x2,0);
end; //function Mahalanobis2
|
|
, включая матлаб и статистику говорят, что это невозможно (по правилу умножения матриц).