Распределение коэффициентов дискриминантной функции

[AndreyL]

Именно так и делается в дискриминантном анализе. А как нужно считать ковариационную матрицу для использования способа, предлагаемого Вами?

[Andrey_Kireew]

AndreyL в сообщении #1427114 писал(а):

Именно так и делается в дискриминантном анализе

Действительно, теперь уже ошибка у меня. Приведённые выше оценки доверительных интервалов неверны, и
по аналогии с регрессионным анализом, как вы предполагали

AndreyL в сообщении #1422313 писал(а):

Есть подозрение, что компоненты вектора связаны с распределением Стьюдента (по аналогии с регрессионным анализом)

получить их не представляется возможным. По той простой причине, что эти методы исходят совсем из разных предпосылок, и распространять одно на другое некорректно.

P.S.: Это вовсе не означает, что теоретические оценки доверительных интервалов получить невозможно. Но подход нужен совсем другой. А с учётом всех "достоинств" дискриминанта Фишера, сама необходимость построения этих интервалов представляется сомнительной.

[AndreyL]

Получил я, похоже, ковариационную матрицу оценок параметров. Если обобщенная ковариационная матрица исходных дискриминаторов $S=\frac{S_x(n_x-1)+S_y(n_y-1)}{n_x+n_y-2}$, то

$$\sigma^2=\frac{\left( n_x+n_y \right)^3-4 \left( n_x+n_y \right) \left( n_x+n_y-1 \right)}{n_x n_y \left( n_x+n_y -k-1 \right) } + \frac{n_x+n_y-2}{n_x+n_y-k-1} \left(\bar{Y}-\bar{X} \right) S^{-1}\left(\bar{Y}-\bar{X} \right)^T$$

при условии $n_x+n_y>k+1$. Тогда, если $W=\left( \begin{array}{cc} 1_{nx} & X \\ 1_{ny} & Y \\ \end{array} \right)$, где $1_{nx}$ и $1_{ny}$ - векторы-столбцы из единичек длиной $n_x$ и $n_y$ соответственно, то ковариационная матрица оценок параметров

$$C=\sigma^2 \left( W^T.W \right)^{-1}$$

[Andrey_Kireew]

AndreyL в сообщении #1440452 писал(а):

Получил я, похоже, ковариационную матрицу оценок параметров.

Поздравляю! идею Вы поняли.
В прошлый раз, я тоже заинтересовался этим вопросом, и получил кое какие формулы. Но, признаться, они отличаются от Ваших, хотя у меня всё выражалось через центрированные факторы, может поэтому. По крайней море, константный член в выражении экв. дисперсии там такой же как у Вас

$$\frac{(n_a+n_b)^3}{n_an_b(n_a+n_b-k-1)}$$

В остальном, различия большие.

Самое неприятное, после многочисленных численных экспериментов выяснилось, что эта формула не работает как должна. Вычисляемая по ней дисперсия оказывается оптимистически заниженной. И это в идеальных условиях.

Если Вы свою формулу не проверяли, то советую проверить

[AndreyL]

А можно ли посмотреть Вашу формулу?
Я пытался делать и через квадратичную форму на ковариационной матрице совмещенных данных, т.е. если $Z=\left( \begin{array}{cc} X \\ Y \\ \end{array} \right)$, то ее ковариационная матрица будет $G=\frac{S_x(n_x-1)+S_y(n_y-1)+\frac{n_x n_y}{n_x+n_y}\left(\bar Y-\bar X \right)^T.\left(\bar Y-\bar X \right)}{n_x+n_y-1}$, но тогда зависимость $\sigma^2$ от $\left(\bar Y-\bar X \right) G^{-1} \left(\bar Y-\bar X \right)^T$ почему-то получалась нелинейной.

-- Чт фев 20, 2020 4:45 am --

И еще вопрос - в какой программе Вы проводите численные эксперименты? Я работаю в Математике и в Матлабе.

[Andrey_Kireew]

AndreyL в сообщении #1440508 писал(а):

но тогда зависимость $\sigma^2$ от $\left(\bar Y-\bar X \right) G^{-1} \left(\bar Y-\bar X \right)^T$ почему-то получалась нелинейной

Да, у меня тоже. По идее, все формулы должны сходится, если это не так, то либо где то ошибка, где то недочёт. Но я помнится много времени на это потратил, и всё получалось складно.

Эксперименты в Matlab. Генерировал пары многомерных нормальных распределений с заданными матожиданиями, делил их на равные части, кажется 100 по 1000, и вычислял 100 векторов коэффициентов. Потом находил их ков. матрицу.
Расхождения были всё время, но не очень большие, но устойчивые, в рамки стат. погрешности они не укладывались. Мне тогда показалось, что это не ошибка, а дефект самой модели, возможно его даже можно скорректировать, но как - пока не знаю.

Попробуйте со своей формулой. может она даже лучше, кто его знает, потом сравним.
Сейчас у меня под рукой их полностью нет, прошло уже 3 месяца, и я как то не думал, что они могут понадобиться. По новой выводить пока нет особого желания, потом, как будет возможность - попробую отыскать.

[AndreyL]

Для начала - зависимость между $\left(\bar Y-\bar X \right) G^{-1} \left(\bar Y-\bar X \right)^T$ и $\left(\bar Y-\bar X \right) S^{-1} \left(\bar Y-\bar X \right)^T$ действительно нелинейная

Используется синтаксис Matlab M

function testSqForms
%% prepare sample sizes and true parameter values
k=4;
nX=5;
nY=7;
meX=2*rand(1,k)-1;
meY=2*rand(1,k)-1;
cov0=cov(2*rand(k+1,k));

%% generate square forms
L=100;
res=ones(L,2);
for i=1:L
 res(i,:)=func(k,nX,nY,meX,meY,cov0);
end;
plot(res(:,1),res(:,2),'.')

function r=func(k,nX,nY,meX,meY,cov0)
%% generate X and Y, calculate main parameter estimation
r=ones(1,2);
X=mvnrnd(meX(ones(nX,1),:),cov0);
Y=mvnrnd(meX(ones(nY,1),:),cov0);
mX=mean(X);
mY=mean(Y);
covX=cov(X);
covY=cov(Y);
S=((nX-1)*covX+(nY-1)*covY)/(nX+nY-2);
r(1,1)=(mY-mX)*inv(S)*(mY-mX)';
G=((nX-1)*covX+(nY-1)*covY+nY*nY/(nX+nY)*(mY-mX)'*(mY-mX))/(nX+nY-1);
r(1,2)=(mY-mX)*inv(G)*(mY-mX)';

Использую цикл for поскольку забыл, как это делается в Матлабе

[Andrey_Kireew]

AndreyL в сообщении #1440518 писал(а):

Для начала - зависимость между ... Использую цикл for поскольку забыл, как это делается в Матлабе

Надо так понимать, что обсуждаемую здесь формулу Вы так и не получили, и экспериментировать Вам пока не с чем. Ну так получайте, ведь как бы там ни было, правилами форума запрещается выкладывать готовые решения.

[AndreyL]

Andrey_Kireew в сообщении #1440539 писал(а):

Надо так понимать, что обсуждаемую здесь формулу Вы так и не получили, и экспериментировать Вам пока не с чем. Ну так получайте, ведь как бы там ни было, правилами форума запрещается выкладывать готовые решения.

Странно, формулу для получения ковариационной матрицы оценок параметров дискриминантной функции я уже выписал в явном виде - мне выкладывать готовые решения не запрещается, поскольку я задаю вопрос, а не отвечаю, и ответ показываю на предмет проверки - правильно/неправильно.
Ну да ладно. Вот скрипт с Монте-Карло. Больше удивляет даже не то, что что дисперсии не сходятся, а то, что оценки параметров оказываются смещенными.

Используется синтаксис Matlab M

function testParameterEstimation
%% prepare sample sizes and true parameter values
k=4;
nX=5;
nY=7;
meX=2*rand(1,k)-1;
meY=2*rand(1,k)-1;
cov0=cov(2*rand(k+1,k));

%% calculate true discriminant function parameters
V=ones(1,k+1);
V(1,2:k+1)=(meY-meX)/cov0;
V(1,1)=V(1,2:k+1)*(meY+meX)'/2;

%disp(func(k,nX,nY,meX,meY,cov0,V));

%% Monte-Carlo
L=500;
res=ones(L,k+1);
for i=1:L
 res(i,:)=func(k,nX,nY,meX,meY,cov0,V);
end;
probab=((1:L)-0.5)/L;

%% plot results
colors=['b.';'c.';'g.';'m.';'y.';'bo';'co';'go';'mo';'yo'];
for j=1:k+1 % number of estimated parameter
 plot(sort(res(:,j)),probab',colors(j,:),'DisplayName',int2str(j));
 hold on;
end;
P = tcdf(sort(res(:,j)),nX+nY-k-1);
plot(sort(res(:,j)),P,'r-','DisplayName','stud');
legend('show')
hold off;

function r=func(k,nX,nY,meX,meY,cov0,param0)
% generate X and Y, calculate main parameter estimation
r=ones(1,k+1);
X=mvnrnd(meX(ones(nX,1),:),cov0);
Y=mvnrnd(meY(ones(nY,1),:),cov0);
mX=mean(X);
mY=mean(Y);
covX=cov(X);
covY=cov(Y);
S=((nX-1)*covX+(nY-1)*covY)/(nX+nY-2);

% calculate discriminant function parameters estimations
Vc=ones(1,k+1);
Vc(1,2:k+1)=(mY-mX)/S;
Vc(1,1)=Vc(1,2:k+1)*(mY+mX)'/2;

% calculate covariance matrix
sigma=((nX+nY)^3+4*(nX+nY)*(nX+nY-1))/(nX*nY*(nX+nY-k-1));
sigma=sigma+(nX+nY-2)/(nX+nY-k-1)*(mY-mX)/S*(mY-mX)';
W=[ones(nX,1),X;ones(nY,1),Y];
C=sigma*inv(W'*W);

% calculate Student stasistics for discriminant function parameters estimations
std=sqrt(diag(C))';
r=(Vc-param0)./std;

Да, кстати, все ли понятно в скрипте, или как-то изменить комментарии?

[Andrey_Kireew]

AndreyL в сообщении #1440562 писал(а):

Странно, формулу для получения ковариационной матрицы оценок параметров дискриминантной функции я уже выписал в явном виде - мне выкладывать готовые решения не запрещается, поскольку я задаю вопрос, а не отвечаю, и ответ показываю на предмет проверки - правильно/неправильно

Возможно Вы меня не поняли, но это было написано в ответ на Вашу просьбу показать мою формулу. Откуда мне знать, получили вы её на самом деле или нет? Если нет, то выложив своё решение, я избавлю Вас от необходимости самостоятельного творчества. Именно это, по идее, форум делать и запрещает.

AndreyL в сообщении #1440562 писал(а):

Больше удивляет даже не то, что что дисперсии не сходятся, а то, что оценки параметров оказываются смещенными.

в Вашем скрипте заложена та же самая формула, которую Вы представили выше. Возможно она и не правильная.

В вашем скрипте, как я понял, на основании заданных средних и ков. матрицы, сначала оцениваются "истинные" параметры, затем генерируется выборка, и по ней однократно оцениваются те же самые параметры. Я бы посоветовал генерировать выборку многократно, и тогда Вы увидите ни какую то эфемерную t - статистику, а реальное распределение параметров. По нему можно вычислить и их дисперсию. Её и сравнивать с теоретической оценкой (она должна мало меняться от выборке к выборке). Это лучше, чем сравнивать единичные значения и потом гадать на кофейной гуще.
В общем, эксперимент Ваш мне не нравится, какой то он премудрый, возможно даже и некорректно поставлен.

[AndreyL]

Andrey_Kireew в сообщении #1440574 писал(а):

В вашем скрипте, как я понял, на основании заданных средних и ков. матрицы, сначала оцениваются "истинные" параметры, затем генерируется выборка, и по ней однократно оцениваются те же самые параметры. Я бы посоветовал генерировать выборку многократно, и тогда Вы увидите ни какую то эфемерную t - статистику, а реальное распределение параметров. По нему можно вычислить и их дисперсию. Её и сравнивать с теоретической оценкой (она должна мало меняться от выборке к выборке). Это лучше, чем сравнивать единичные значения и потом гадать на кофейной гуще.
В общем, эксперимент Ваш мне не нравится, какой то он премудрый, возможно даже и некорректно поставлен.

Теперь я Вас не очень понимаю. Тогда объясню свою тактику: сначала генерируются два вектора истинных средних и истинная ковариационная матрица. Именно для этих параметров и считаются истинные значения параметров дискриминантной функции. Потом 500 раз генерируются две выборки (их объемы всегда одинаковые) случайных величин распределенных нормально с известными центрами и известной одной и той же ковариационной матрицей. Естественно, что оценки центров и ковариационных матриц по выборкам будут отличаться от истинных значений. По этим оценкам считаются оценки параметров дискриминантной функции и их ковариационная матрица (и так 500 раз). Далее для каждого параметра считается t-статистика. Просто генерация выборок и дальнейшая работа с ними убрана в отдельную функцию, которая вызывается 500 раз.
А как Вы предлагаете?

[Andrey_Kireew]

Я Вам советую объединить все сгенерированные коэффициенты в одну матрицу (можно их возвращать из той же функции) и напрямую вычислить ковариационную матрицу коэффициентов. Потом её можно сравнить с усреднённой ковариационной матрицей, вычисленной по Вашей формуле

[Andrey_Kireew]

AndreyL в сообщении #1440562 писал(а):

Больше удивляет даже не то, что что дисперсии не сходятся, а то, что оценки параметров оказываются смещенными

И на каком основании Вы делаете такие выводы? С чем не сходятся дисперсии? Друг с другом? И как Вы это проверили, если Они вычисляются внутри функции, и там же остаются? по t критерию?

Ну и смещение коэффициентов как Вы выявили? они же у Вас тоже не выходят за пределы функции.

Может, конечно, у Вас для этого есть другой скрипт, но тогда не совсем понятно, зачем Вы привели этот, от которого всё равно нет никакого толку

[AndreyL]

Andrey_Kireew в сообщении #1440587 писал(а):

Я Вам советую объединить все сгенерированные коэффициенты в одну матрицу (можно их возвращать из той же функции) и напрямую вычислить ковариационную матрицу коэффициентов. Потом её можно сравнить с усреднённой ковариационной матрицей, вычисленной по Вашей формуле

Не совсем понимаю, как это можно сделать. При каждой реализации случайных выборок $\bar Y-\bar X$ и $S$ разные, поскольку они по этим выборкам и считаются. Тогда и дисперсия оценки в каждой реализации своя. Или лучше использовать истинные дисперсии?

[Andrey_Kireew]

AndreyL в сообщении #1440799 писал(а):

При каждой реализации случайных выборок $\bar Y-\bar X$ и $S$ разные, поскольку они по этим выборкам и считаются. Тогда и дисперсия оценки в каждой реализации своя

Да, оценка дисперсии по формуле каждый раз разная, а по набору коэффициентов Вы получите только одну ков. матрицу. Но это нормально. Постройте гистограмму распределения теоретических оценок и отметьте на ней же дисперсию, полученную непосредственно по коэффициентам. Если формула корректная, то оценка дисперсии, сделанная по коэффициентам должна соответствовать наиболее вероятному значению гистограммы, и быть в пределах выборочной ошибки дисперсии (как её найти для своей выборки, думаю, найдёте сами). Заодно проверите, соответствует ли полученное распределение распределению Хи-квадрат.
При этом, можно ещё построить распределения коэффициентов, и посмотреть, насколько оно соответствует нормальному.

Насколько я понимаю, Ваша цель, определить насколько хорошо согласуется оценка дисперсии, полученная по формуле, с реальной дисперсией. По t критерию можно только узнать, соответствует ли распределение res распределению Стьюдента с соответствующим числом степеней свободы, или нет. Сравнивать нужно тоже по Хи-квадрат. И в итого, всё что вы получите - это согласуются дисперсии оценок, или нет. Насколько они различаются, и смещены ли они - вообще сказать ничего нельзя. А Вы ещё зачем то начинаете вычислять по res вероятности. Задайте себе вопрос, вероятности чего Вы вычисляете? или вам это не важно. А в итоге строите распределение этих вероятностей. Что Вы надеетесь увидеть по этому распределению? Лично мне всё это представляется как бессмысленное нагромождение стат. методов. Обычно так часто делают, чтобы произвести впечатление, и придать наукообразия своим работам. Но здесь то все участники анонимные, в том числе и Вы сами, какой в этом смысл?