Open kel-github opened 4 weeks ago
crnolan
commented
2 weeks ago Unless I've misunderstood something, BrentZero appears to be a method to do better than a bisection search algorithm to find the zero point of a function. I would assume that we can instead use an existing search algorithm already implemented in R and only resort to replicating BrentZero if we are having performance issues.
crnolan
commented
2 weeks ago ... which might not even matter, as the uniroot function already finds the zeros in a specified range for a defined function, and references Brent (1973).
kel-github
commented
2 weeks ago List of sensible order for implementation of functions:
jorittmo
commented
2 weeks ago Just as an addition, the p-value given from Roy's statistic as implemented in base (stats) R manova() is using the f-distribution and only gets an upper bound
The most magical thing about Psy is that it has functions to compute the largest GCR of a Wishart matrix that you could expect to obtain by chance, given the size of your covariance matrix. We need to work out how best to implement these functions for an R package.
What PSY does to compute the largest GCR that can be expected under the null hypothesis:
Note: the below paper resources can be found in the folder computing-largest-GCR-under-the-null which is in the zotero library here.
All pascal files and associated comments are in this folder of the repository.
The critical value for the GCR, i.e. the largest GCR (theta) that could be expected under the null hypothesis) is broadly computed as follows:
Step 1. compute p(theta <= x) using a numerical approximation method (see step 1 below)
Step 2. evaluate the p(theta <= x) across a range of x's, until you find the x where p(theta <= x) == .975. This X is the critical value for the GCR.
Step 1: for any given theta (root value), compute the probability that theta is less than or equal to some given value (x), given the number of variables - aka the size of the covariance matrix - i.e. p(theta <= x)
The goal in this step is to compute eq 3 of Pillai (1965). On the distribution of the largest characteristic root of a matrix in multivariate analysis. This equation states how to compute p(theta <= x).
See 'GCR.pas' and 'PsyBaseStats.pas' for the relevant functions (listed below). See Bird (2002) in the 'logic-of-...' folder of the zotero library for the definition of parameters s, m, & n.
The key functions to achieve this are: Ksmn: this computes the constant factor in the largest root CDF. See Pillai (1965) On the distribution of the largest characteristic root of a matrix in multivariate analysis. Eq 2.
RoyExact: This is an exact expression for obtaining p(theta <= x) - See Pillai (1954) On some distribution problems in multivariate analysis.. Eqs 2.3.2, 2.3.3 & 2.3.4.
RoyExact calls BetaZero and BetaRoy from PsyBaseStats.pas, which computes the beta functions defined in Pillai (1954). Also see Pillai (1965) Eq 5.
BetaRoy calls the IncBetaContFrac which is the meat of the implementation of the beta function as mentioned above. This is based on the numerical approximation method for computing the first and second derivatives of an incomplete beta function outlined in Boik & Robinson-Cox (1999). Derivatives of the incomplete beta function. This paper shows that the incomplete beta function can be written as a hypergeometric series. Thus, computing eqs 2.3.2, 2.3.3, & 2.3.4 in Pillai (1954) requires summing over the relevant series as defined by Boik & Robinson-Cox (1999).
PillaiApproxOriginal allows one to approximate the p(theta <= x) when the value of the s parameter does not allow use of the exact solution. This function is used when s is between 4 and 20. See Pillai & Flury (1984). Percentage points of the largest characteristic root of the multivariate beta matrix Eq 2.3. Note that the approximation method is improved (i.e. attenuating the impact of misrepresentations in floating point values) by using the double sum method. See function CompSum and wikipedia.
RoyExact and PilaiApproxOriginal are called by gcr_CritProb_Exact and gcr_CritProb_Approx which performs the subtraction between p(theta <= x) and the desired probability (e.g. .05). These latter two functions are called by gcr_Crit which assigned the appropriate function for which zero needs to be found (see below).
Step 2: evaluate p(theta <= x) across a range of x values, until you find p(theta <= x) == .05.
gcr_Crit calls BrentZero which takes as its input a function and the associated arguments. It then finds the value of x (in this case, the critical GCR value we want to find) that returns zero as the output when the p(theta <= x) - critical p (e.g. .05) = 0. See Brent (1971). An algorithm with guaranteed convergence for finding a zero of a function