Building PDFs of an eigenvalue of Wishart distribution in MATLAB

7 次查看（过去 30 天）
JoonYoung
JoonYoung 2022-9-20
评论： Bjorn Gustavsson 2022-9-22
I have been trying to build MATLAB code from PDFs in a paper, "On the marginal distribution of the eigenvalues of Wishart matrices." Specifically, I am working on
This is the PDF of the largest eigenvalue of the uncorrelated central Wishart matrix
This is the PDF of a generic (unordered) eigenvalue of the uncorrelated central Wishart matrix
which are Equation (38) in the paper. And I have constructed a MATLAB script for these two PDFs like this.
input_end=30; % Define the length of time domain
x=0:0.01:input_end; % Define time domain
y_lmda1=zeros(1,length(x)); % Define f_lmda1(x)
y_generic=zeros(1,length(x)); % Define f_generic(x)
M=5;
N=5;
nmax=max(M,N);
nmin=min(M,N);
%%% Define K %%%
product_of_K_when_max=1;
for i=1:nmin
fact_of_prdct_nmax=factorial(nmax-i);
product_of_K_when_max=product_of_K_when_max*fact_of_prdct_nmax;
end
product_of_K_when_min=1;
for i=1:nmin
fact_of_prdct_nmin=factorial(nmin-i);
product_of_K_when_min=product_of_K_when_min*fact_of_prdct_nmin;
end
K=(product_of_K_when_min.*product_of_K_when_max)^(-1);
%%% Define the PDF of the largest eigenvalue %%%
PDF_sum=0; %%% Define the initial state of the matrix omega_lambda1 %%%
PDF_generic_doublesummation=0;
omega_lambda1=zeros(nmin,nmin); %%% Allocating initial values to the matrix Omega%%
%%% Define the matrix omega_lambda1 %%%
omega_lambda1=zeros(nmin,nmin); %%% Allocating initial values to the matrix Omega %%%
for k=1:length(x); %% output will be calculated with corresponding x(k)%%
for n=1:nmin;
for m=1:nmin;
for i=1:nmin;
for j=1:nmin;
if i<n && j<m
omega_lambda1(i,j)=gamma(i+j-2+nmax-nmin+1)*gammainc(x(k),i+j-2+nmax-nmin+1);
%%% gamma() is multiplied because the incomplete Gamma function in MATLAB and the incomplete Gamma function in the paper are different.%%%
elseif i>=n && j >= m
omega_lambda1(i,j)=gamma(i+j+nmax-nmin+1)*gammainc(x(k),i+j+nmax-nmin+1);
else
omega_lambda1(i,j)=gamma(i+j-1+nmax-nmin+1)*gammainc(x(k),i+j-1+nmax-nmin+1);
end
end
end
PDF_sum=PDF_sum+((-1).^(n+m)).*(x(k).^(n+m-2+nmax-nmin)).*(exp(-x(k))).*(det(omega_lambda1));
omega_lambda1=zeros(nmin,nmin);
end
end
y_lmda1(k)=K.*PDF_sum; %%% After double summations, the normalizing factor is multiplied and it's input to the output vector.%%%
PDF_sum=0; %%% Initializing the summation loop with new x value. %%%
end
%%% Define the matrix omega_generic %%%
sigma_nmin_1=nmin.^(-1/(nmin-1));
omega_generic=zeros(nmin,nmin);
for k=1:length(x);
for n=1:nmin;
for m=1:nmin;
for i=1:nmin;
for j=1:nmin;
if i<n && j<m
omega_generic(i,j)=factorial(i+j-2+nmax-nmin)*sigma_nmin_1;
elseif i>=n && j >= m
omega_generic(i,j)=factorial(i+j+nmax-nmin)*sigma_nmin_1;
else
omega_generic(i,j)=factorial(i+j-1+nmax-nmin)*sigma_nmin_1;
end
end
end
PDF_generic_doublesummation=PDF_generic_doublesummation+...
((-1).^(n+m)).*(x(k).^(n+m-2+nmax-nmin)).*(exp(-x(k)))*(det(omega_generic));
omega_generic=zeros(nmin,nmin);
end
end
y_generic(k)=K.*PDF_generic_doublesummation; % After double summations, the normalizing factor is multiplied and it's input to the output vector.
PDF_generic_doublesummation=0; %%% Initializing the summation loop with new x value.%%
end
figure;
plot(x,y_lmda1,'k-.',x, y_generic,'k--');
legend("lambda1", "generic",'FontSize', 14);
ylabel('f(x)');
xlabel('x');
%%% PDF verification %%%
sum(y_lmda1)
sum(y_generic)
I normalized the PDFs with K. However, my result is here , which is not valid because both PDF values go way over 1. I have been trying to figure out what is the problem but it seems it is beyond my capability. I want to know what I am missing now from others.
  1 个评论
Bjorn Gustavsson
Bjorn Gustavsson 2022-9-22
Since the distributions seems to be continuous shouldn't the normalization be done by integration rather than summation?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by