function [t s] = TW_trace_ratio_threshold(p,n,beta,alpha)
%function [t s] = TW_trace_ratio_threshold(p,n,beta,alpha)
% Let L_j be the eigenvalues of the sample covariance matrix of multivariate
% Gaussian obserations with n samples in p dimensions,
%
% INPUT: p = dimension, n = number of samples,
% beta = 1/2 real/complex valued data,
% alpha = desired false alarm / right tail probability
%
% OUTPUT: t = An approximate value such that Pr[ L_1/(1/p sum(L_j) )> t ] ~ alpha
% s = Pr[ (L_1/(1/p sum(L_j) - mu)/sigma > s] ~ alpha
%
%
% Threshold s is found by interpolation of precomumputed values of the
% Tracy-Widom distribution, and its second derivative.
% These were generated using code by Prof. Folkmar Bornemann, based on his paper
% On the Numerical Evaluation of Distributions in Random Matrix Theory, 2010.
%
if beta==1
load TW_beta1.mat
else
load TW_beta2.mat
end
% mu_np and sigma_np - centering and scaling
[mu_np sigma_np] = KN_mu_sigma(n,p,beta); mu_np = mu_np ./ n; sigma_np = sigma_np ./ n;
% Eq. 1.7 in paper [1]
U_cdf_complementary = 1-TW_s +1/beta/n/p * (mu_np/sigma_np)^2 * TW_s_tag_tag;
% next find index where probability is close to alpha
[val idx] = min(abs(U_cdf_complementary-alpha) );
t = mu_np + sigma_np * x(idx);
s = x(idx);