Abstract:
A method of classifying high-dimensional data to facilitate subsequent processing of the data. Generally characterized, the method includes determining a statistical distribution function for a first dataset; determining a statistical distribution function for a second dataset; identifying a characterizing feature of the first and second datasets; determining a probability of distinction based on the first and second statistical distribution functions, with the probability of distinction being the probability that the characterizing feature can be used to distinguish the first dataset from the second dataset; and identifying a subject of the first dataset based on the probability of distinction. Once the subject has been identified, the first dataset can be further analyzed using processing techniques appropriate to the subject.

Description:
RELATED APPLICATIONS  
       [0001]     The present U.S. non-provisional patent application claims priority benefit of an earlier-filed provisional patent application titled METHODS FOR ANALYZING DATA, Ser. No. 60/802,365, filed May 22, 2006. The identified earlier-filed application is hereby incorporated by reference into the present application. 
     
    
     FIELD OF THE INVENTION  
       [0002]     The present invention relates generally to biostatistical methods of classifying data, and, more particularly, to a method of classifying high-dimensional data to facilitate subsequent processing of the data.  
       BACKGROUND OF THE INVENTION  
       [0003]     Cancer, the abnormal and uncontrolled division of cells causing benign or malignant growth, is the second leading cause of death in the U.S., exceeded only by heart disease. The early detection, classification, and prognosis of cancer remains the best way to ensure patient survival and quality of life.  
         [0004]     Many current techniques to predict and classify cancer use deoxyribonucleic acid (DNA) microarray data, and, in particular, gene expression data, which allows for testing large quantities of genetic material. Generally, using the data, cells can be analyzed, genetic information can be extracted, and phenotypes can be determined, i.e., the characteristics of an organism, as determined by its genes and the relationship of the genes to the environment.  
         [0005]     Typically, DNA molecules are placed on a microarray and attached using probes. Each of these DNA molecules may include a complement, generally referred to as a target or messenger ribonucleic acid (mRNA), which is, generally speaking, the working copy of genes within cells. When testing a cell from an organism, the mRNA is added to the microarray where the mRNA “hybridizes” with the DNA to which it is a complement. Thus, some of the discrete spots will contain mRNA hybridized to DNA, while other spots will not. In some technologies, the cellular mRNA is reverse-transcribed onto complementary DNA (cDNA), which is a complementary copy of the mRNA. The cDNA is linearly amplified and is subsequently used to hybridize with the probes.  
         [0006]     The marker attached to the mRNA is used to determine which spots contain mRNA hybridized to DNA. Typically, the markers are fluorescent molecules that fluoresce when exposed to laser light of an appropriate frequency and power. This fluorescence can be measured to determine the degree to which a gene has expressed. For example, a high fluorescence detected for a particular gene indicates that the gene is very active. Conversely, a low or no fluorescence for a gene indicates that the gene is inactive. Thus, by examining the fluorescence of the microarray, the degree of activity of the different genes can be determined.  
         [0007]     Using current techniques, it is possible to determine gene function for tissues affected and unaffected by a disease. By comparing the two phenotypes, it can be determined which genes contribute to certain diseases and how they contribute. However, current techniques are limited. In particular, the detection and classification of gene expression data is based on only a few observations. In many cases, only a few underlying gene components account for much of the data variation. For example, only a few linear combinations of a subset of genes may account for nearly all of the response variation. Unfortunately, it is exceedingly difficult to determine which genes are members of the subset given the number of genes in the microarray and the small number of observations. As such, given the number of genes on a single microarray (usually is in the thousands), a relatively small number of observations cannot provide accurate statistical data.  
       SUMMARY OF THE INVENTION  
       [0008]     The present invention provides a method of classifying high-dimensional data to facilitate subsequent processing of the data. In one embodiment, the method comprises the steps of determining a statistical, e.g., cumulative, distribution function for a first dataset; determining a statistical distribution function for a second dataset; identifying a characterizing feature of the first and second datasets; determining a probability of distinction based on the first and second statistical distribution functions, wherein the probability of distinction is the probability that the characterizing feature can be used to distinguish the first dataset from the second dataset; and identifying a subject of the first dataset based on the probability of distinction. Once the subject has been identified, the first dataset can be further analyzed using processing techniques appropriate to the subject.  
         [0009]     The first and second datasets may include microarray data, in which case the characterizing feature may be an expression of a gene associated with a particular biological mechanism such as, e.g., cancer. Alternatively, the first and second datasets may include surface-enhanced laser desorption ionization-time of flight data, in which case the characterizing feature may be a mass-to-charge ratio. The method may further include the step of filtering the data to eliminate potential characterizing features which would have low probabilities of distinction. The filter used in filtering the data may be based at least in part on the nature of the data or on the recognition of patterns in the data.  
         [0010]     The method may further include the steps of identifying a plurality of characterizing features; determining a probability of distinction for each of the characterizing features; and selecting whichever of the characterizing features has the highest probability of distinction.  
         [0011]     The method may further include the step of using the probability of distinction to specify a distribution of hyperparameters in a hierarchical model such as, e.g., a Bayesian neural network or a Gaussian process model.  
         [0012]     It will be appreciated by those with ordinary skill in the art that, among other advantages, the present invention allows for, when genes are the subject of the data, determining a probability of importance for each gene. This probability, combined with some basic biological knowledge, may provide a better understanding of the genes, and, therefore, may improve the detection, classification, and diagnosis of diseases such as cancer.  
         [0013]     These and other features of the present invention are described in greater detail in the section titled DETAILED DESCRIPTION, below. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0014]     The following drawings are included to further demonstrate and clarify certain aspects of the present invention and are not meant to be limiting:  
         [0015]      FIG. 1  is a table of features, associated samples, and associated probabilities of distinction;  
         [0016]      FIG. 2  is an exemplary algorithm of a bootstrap technique;  
         [0017]      FIG. 3  is a scatter plot of probabilities of selecting an m/z ratio interval;  
         [0018]      FIG. 4  is a table of bootstrapped p-values for SFS of m/z intervals;  
         [0019]      FIG. 5  is a comparison between a traditional method and a Gaussian process method;  
         [0020]      FIG. 6  is a plot of an h-surface generated network with a 2-dimensional input, in accordance with embodiments of the disclosure;  
         [0021]      FIG. 7  is a table of analyzed ovarian cancer data;  
         [0022]      FIG. 8  is a table of features at various significance levels and thresholds;  
         [0023]      FIG. 9  is a plot of a simulation of hidden-to-output weights;  
         [0024]      FIG. 10  is a plot of a simulation of five input-to-hidden weights;  
         [0025]      FIG. 11  is a table comparing various classifier methods; and  
         [0026]      FIG. 12  is a listing of source code for implementing a method for analyzing data.  
     
    
     DETAILED DESCRIPTION OF THE INVENTION  
       [0027]     With reference to the accompanying figures, a method is herein described, shown, and otherwise disclosed in accordance with one or more preferred embodiments of the present invention.  
         [0028]     In one respect, the present invention provides methods for classifying high-dimensional data, and in particular, microarray data or surface-enhanced laser desorption ionization-time of flights (SELDI-TOF) data. The methods may use a “shallow feature selection” (SFS) that assigns each feature a probability of being selected based on, for example, an empirical distribution of the data. A “feature”, as used herein, refers to an attribute or characteristic of the data. For example, for data concerning tissue, a feature of the data may be a protein expression level of the tissue. A “probability of distinction”, or simply a “probability”, as used herein, refers to a probability of a selected feature being able to distinguish between data classes (e.g., a normal sample and a cancerous sample). Using these methods, the microarray data or SELDI-TOF data can be reduced to an applicable case for subsequent processing.  
         [0029]     In another respect, the probability assigned to a given feature may be used to specify a distribution or distributions of hyperparameters in hierarchical models of, for example, Bayesian neural networks (BNN) or Gaussian process models. Neural networks, such as BNN models, can be used in the classification and regression of data due to their ability to model complex non-linear systems. Theoretically, neural networks can simulate any continuous function in a compact range. Neural networks can simulate a continuous function if the number of hidden units approaches infinity, where “hidden units”, as used herein, refers to the connection units between an input layer and an output layer.  
       Shallow Feature Selection  
       [0030]     A normalized training set, D, with n samples is equal to
 
{( x   u   , t   u )| x   u ε           m   , t   u ε{0,1} and μ=1, 2, . . .  n}   Eq. 1
 
 where t u ε{0,1} are class labels and samples x u  are recorded in column vectors. In some embodiments, feature selection may be based on pattern recognition. Additionally or alternatively, some characteristics of the dataset itself may be used to filter useless features, such as features that have the same distributions across two different classes. 
 
         [0031]     For example, referring to  FIG. 1 , each sample, x i , is denoted by a column vector, where f i  represents a feature or attribute, e.g., the expression of a gene in microarray data or the mass-charge ratio in surface-enhanced laser desorption ionization-time of flight (SELDI-TOF) data, and p i  is the probability a feature is selected. f i  may be useful if the empirical distributions of the two datasets are statistically distinct. In other words, f i  may be useful for two datasets with different data distributions. Moreover, for a reasonable representation of the training data, each feature in question may be associated with a prior probability of being selected for classification, independent of particular classifiers. These features, as used herein, refer to shallow feature selections (SFS).  
         [0032]     Without any assumption of population distribution, the empirical distributions of values in X i =(x i,1 , . . . , x i,k ) and X′ i =(x i,k+1 , . . . , x i,n ) may be compared using, for example, a two-sided Kolmogorov-Smirnov (KS) test (i.e., the null hypothesis H o  where X i  and X′ i  are drawn from the same continuous distribution) with a given level of significance, α.  
         [0033]     Due to the high dimensionality and relatively small amount of training data, the statistical difference between the control dataset {x 1 , . . . , x c } and the treatment dataset {x c+1 , . . . x n } may be hard to survey. In one respect, f i  should be assumed to be more important if the control dataset {C i ={x i,1 , . . . , x i,c }) and the training dataset (T i ={x i,c+1 , . . . , x i,n }) are more probably from distinct populations.  
         [0034]     Assuming X 1 , . . . X c   iid ˜F 1 (z) and Y 1 , . . . , Y t   iid ˜F 2 (z), the cumulative distribution function (CDF) F 1 (z) can be approximated by the empirical CDF {circumflex over (F)} 1 (z). By Gnedenko theorem, if F 1 (z)=F 2 (z), the statistic of the D c,t =sup|{circumflex over (F)} 1 (z)−{circumflex over (F)} 2 (z)| tends to be 0 as c,t approaches infinity, ∞.  
         [0035]     Moreover, Kolmogorov-Smironov theorem affirms that  
             ct   n       ⁢     D     c   ,   t         ∼     K   ⁡     (   z   )           
 
 where  
         K   ⁡     (   z   )       =     1   -     2   ⁢       ∑     i   =   1     ∞     ⁢           ⁢         (     -   1     )       i   -   1       ⁢     exp   ⁡     (       -   2     ⁢     i   2     ⁢     z   2       )                   
 
 is the CDF of Kolmogorov&#39;s distribution. Therefore, asymptotically, D c,t  may serve as a test statistic for the null hypothesis H o :F 1 (z)=F 2 (z). 
 
         [0036]     The strength of evidence against the null hypothesis H o  is measured by a p-value defined by Pr(T&gt;t obs |H o ) where t obs  is the value of a test statistic, T, for the observed data. In one respect, a smaller p-value is regarded as stronger evidence against H o .  
         [0037]     Referring to  FIG. 2 , an example algorithm based on the bootstrap technique of the present invention assigns a probability to f i  based on the control dataset, C i , and treatment dataset, T i , denoted by Pr(f i |a=a o ).  
         [0038]     The bootstrap method requires no theoretical calculations no matter how mathematically complicated the estimator, {circumflex over (θ)}. Pr(f 1 |a 0 ) may be viewed as the degree of feature f i  distinguishing the two datasets (e.g., control and training) at a preciseness level of α o . In some embodiments, the SFS of the present invention provides users some prior knowledge of the importance of features for classification in view of the data structure itself. Furthermore, the empirical distribution of bootstrapping p-values {circumflex over (θ)} are able to show some statistical characteristics of the considered feature; for example, a small coefficient of variation is preferred.  
         [0039]     After repeating the procedure (N−1) more times on the resampled training datasets, the simple voting will result in Pr(f i |a)=#(H o rejected)/N, the probability of choosing f i  as a distinguishing feature at a level of confidence, α.  
         [0040]     Referring to  FIG. 3 , a graph is shown illustrating the probability of selecting a certain feature. In this example, an ovarian dataset from the National Cancer Institute (NCI) was analyzed. For a particular mass-charge (m/z) ratio, a probability of selecting that ratio is determined. It is noted that while the NCI&#39;s dataset was used, in particular for the relevance of the data to cancer classification, other datasets may also be analyzed.  
         [0041]     In some embodiments, feature selection may be based on the p-values guaranteed by n i n i /(n i +n′ i )≧4, where n i , n′ i  are the sample sizes of X i  and X′ i  respectively. Additionally, the bias and standard error calculated from the N bootstrap (a repeating step) replicates of the p-value may provide more information about SFS (see  FIG. 4 ), especially for the amphibolous case of Pr(f i |)≈0.05, a standard confidence test level generally used in the art. Although the SFS is independent of specific classifiers, it may not isolate from classification. In BNN models, Pr(f i |) may be used to specify the covariance matrix of Gaussian processes and the prior distributions of input-to-hidden weights.  
         [0042]     In contrast to some other classifier-independent approaches to data reduction, such as Principle Component Analysis (PCA), wavelet analysis, etc., SFS is capable of selecting a set of features without losing their biological meaning. In some embodiments, features with large discriminant power to distinguish two different datasets (e.g., a cancer dataset versus a normal dataset) may have a higher likelihood of being selected compared to features with small discriminant power. Furthermore, the probability assigned to each feature is actually a kind of prior knowledge about its contribution to classification, hence use of the term “prior probability” in this disclosure. In contrast, neither PCA nor wavelet utilizes class information in data reduction.  
       Binary Classifications of BNN  
       [0043]     In contrast to conventional training methods for multilayer perceptron (MLP) that are based on maximum likelihood estimation, the Bayesian approach to neural networks utilizes the prior distribution of network weights before training and revises it after observing the data. Therefore, BNN can provide not only the means of predictive weights but also their uncertainties. Another advantage of BNN is the appropriate choices of a number of hidden layers and their dimensions, which are almost impossible to solve by traditional neural networks.  
         [0044]     The inference of an underlying Bayesian binary classifier from D parameterized by the weight vector w, is to find f(x,w|D,M) for the probability of x being in class 1 in model M. In some embodiments, M may be a logistic model where f(x,w)=[1+exp{−h(x,w)}] −1 , satisfying ∫f(x,w)dx&lt;∞ and where h(x, w) is  
               h   ⁡     (     x   ,   w     )       =     log   ⁢       Pr   ⁡     (       t   =     1   ❘   x       ,   w     )         Pr   ⁡     (       t   =     0   ❘   x       ,   w     )                   Eq   .           ⁢   2             
 
         [0045]     Considering the approximation ability, h, is adopted to be the adaptive basis function of one hidden layer (with u units) perceptron h(x,w)=b+W T ψ(b+M T x). Since the biases b and b, if necessary, can be absorbed in M T  and W T , respectively, h(x, w) may be considered as
 
 h ( x,w )= W   T ψ( M   T   x )= W   T   y   Eq. 3
 
 where w=(M, W) and ψ is a bounded activation function, such as hyperbolic tangent tanh(z)=[exp(−z)]/[exp(z)+exp(−z)], or error function  
         erf   ⁡     (   z   )       =       2     π       ⁢       ∫   0   z     ⁢       exp   ⁡     (     -     t   2       )       ⁢           ⁢       ⅆ   t     .               
 
 The binary classification is a noise-free regression problem of logit (Pr(xεclass 1))=h(x, w), where logit(p)=log(p/1−p)). As such, the likelihood is  
               Pr   ⁡     (     D   ❘   w     )       =       ∏     μ   =   1     n     ⁢           ⁢         f   ⁡     (       x   μ     ,   w     )         t   μ       ⁢       (     1   -     f   ⁡     (       x   μ     ,   w     )         )       1   -     t   μ                     Eq   .           ⁢   4             
 
         [0046]     Once given p(w), the posterior distribution of parameters p(w|D) (proportional to p(w)p(D|w) may be clear, and the probability of testing sample x being in class 1 is  
                     Pr   ⁡     (         t     n   ⁢           ⁢   ɛ   ⁢           ⁢   w       =     1   ❘     x     n   ⁢           ⁢   ɛ   ⁢           ⁢   w           ,   D     )       =       ⁢     ∫       f   ⁡     (       x     n   ⁢           ⁢   ɛ   ⁢           ⁢   w       ,   w     )       ⁢     p   ⁡     (     w   ❘   D     )       ⁢     ⅆ   w                     ≈       ⁢       1   N     ⁢       ∑     i   =   1     N     ⁢           ⁢     f   ⁡     (       x     n   ⁢           ⁢   ɛ   ⁢           ⁢   w       ,     w   i       )                         Eq   .           ⁢   5             
 
 where w i  are drawn from a population of p(w|D). The intuitive assumption of p(w)−∝exp{−ξw T w/2} relies on the belief that smaller weights lead to less certain predictions, where the hyperparameter ξ is defined in a hierarchical model, for example, the automatic relevance determination. The integral in Eq. 5 can be evaluated using, for example, a Metropolis algorithm. 
 
         [0047]     Again, the desired p(w) may be viewed as the prior distribution over all the functions described by Eq. 3. Motivated by the Gaussian process (GP) approach to nonparametric Bayesian regression, a large class of neural network models converge to Gaussian processes in the limit of an infinite number of hidden units, and the overfitting problem may be alleviated by some hierarchical model of hyperparameters. Theoretically, finite-dimensional distributions of GP may be determined by its finite mean and covariance functions.  
         [0048]     An advantage of using a GP approach is that there is no need to make complicated assumptions on the distribution of w. For a covariance function, however, a distribution of h(x nεw , w|D) may be achieved. Referring to  FIG. 5 , a comparison between a traditional method and the present GP method is shown.  
         [0049]     Although the logistic model, in terms of latent values h=h(x, w), may be specified by a wide variety of covariance functions, mainly it is in the form of
 
Cov( h,h′ )=δ hh     t     ε+c   2   +x   T  diag(γ 1   2 , . . . γ m   2 ) x′+β   2  exp{−( x−x′ ) T  diag(λ 1   2 , . . . , λ m   2 )( x−x′ )}  Eq. 6
 
 where c 2 , γ i   2 , β 2 , λ i   2  are hyperparameters. Due to the computational issue, the covariance function often adds a “jitter” of δ hh′ ε where ε&gt;0 and δ hh′  is known as the Kronecker delta. The term of c 2 +x T  diag(γ 1   2 , . . . λ m   2 )x′ is for the linear regression, where the last term is based on the assumption that nearby inputs will have highly correlated outputs, allowing a different distance measure for each input dimension than the smaller λ 1   2 , which is less important of i th  features. 
 
         [0050]     The joint distribution of GP h n+1 ,=(h 1 , . . . , h n , h nεw ) T  given X n+1 =(x 1 , . . . x n , x nεw ) and the parameterized covariance function matrix K n+i  is
 
 h   n+1   |X   n+1   , K   n+1   ˜N   n+1 ( Eh   n+1   , K   n+1 )  Eq. 7
 
 where  
           K     n   +   1       =     (             K   n               k   new   T           ⁢           K   new             k           )       ,     
     ⁢       K   n     =       (     Cov   ⁡     (       h   μ     ,     h   v       )       )     mxn       ,     
     ⁢       k   new     =       (       Cov   ⁡     (       h   1     ,     h   new       )       ,   …   ⁢           ,     Cov   ⁡     (       h   n     ,     h   new       )         )     T           
     and     
       k   =         Cov   ⁡     (       h   new     ,     h   new       )       T     .         
 
 With the assumption of Eh=0 for the equal prior probability of each class, the conditional distribution of h new  with hyperparameters gives t n =(t 1 , . . . , t n ) or equivalently h n =(h 1 , . . . , h n ) T , and
 
 h   new   |x   new   , D, K   n+1   ˜N ( k   new   T   K   n   −1   h   n   k k−k   new   T   K   n   −1   k   new   Eq. 8
 
         [0051]     For any hyperparameter in θ in K n , it can be replaced by the maximum likelihood estimate (MLE) from          L/         θ=0, where L is the log-likelihood of h n  and is expressed as  
             L   =         -     n   2       ⁢     log   ⁡     (     2   ⁢   π     )         -       1   2     ⁢   log   ⁢           ⁢   det   ⁢           ⁢     K   n       -       1   2     ⁢     h   n   t     ⁢     K   n     -   1       ⁢     h   n                 Eq   .           ⁢   9                   ϑ   ⁢           ⁢   L       ϑ   ⁢           ⁢   θ       =         -     1   2       ⁢     tr   ⁡     (       K   n     -   1       ⁢       ϑ   ⁢           ⁢     K   n         ϑ   ⁢           ⁢   θ         )         +       1   2     ⁢     h   n   T     ⁢     K   n     -   1       ⁢       ϑ   ⁢           ⁢     K   n         ϑ   ⁢           ⁢   θ       ⁢     K   n     -   1       ⁢     h   n                 Eq   .           ⁢   10             
 
         [0052]     Referring to  FIG. 6 , a simulation of random h-surface for MLP with a 2-dimensional input generated by the MLE-based GP method is shown. The z-axis is the h-surface of a randomly generated network with a 2D input (x 1 , x 2 ) T , using Eq. 6, the x-axis is x 1  and the y-axis is x 2 .  
         [0053]     Referring again to Eq. 6, letting y i  be the i th  component of y=ψ(M T x), and W˜N u (0, σ 0   2 I u ) and Eh=0 the covariance function for Eq. 3 is  
               Cov   ⁡     (     h   ,     h   ′       )       =         σ   0   2     ⁢       ∑     i   =   1     u     ⁢           ⁢     E   ⁡     (       y   i     ⁢       y   ′     i       )           =     u   ⁢           ⁢     σ   0   2     ⁢       E   ⁡     (       y   u     ⁢     y   u   ′       )       .                 Eq   .           ⁢   11             
 
 Assuming that M., 1 , . . . , M., u   iid ˜N m (0,Σ), where Σ=diag(σ 1   2 , . . . , σ m   2 ), the analytic covariance functions for both error and Gaussian transfers coincide with Eq. 6 under some conditions, such as when the variance is larger than the largest distance in x-space. 
 
         [0054]     In some embodiments, parameter tuning may not be necessary for BNN models because of the hidden noisy features in high-dimensional data. Instead of being initialized randomly or trained under strong assumption of distributions, the MLE-based GP model is specified by λ i   2 =γ i   2 =Pr(f i |α), β 2 =c 2 =1, ε=0.001 where Pr(f|a) weakens the influence of “useless” features in the covariance function. This strategy may simplify the GP model and yields an excellent performance on specific gene data, as reflected in the ovarian data analysis shown in  FIG. 7 , where k is the matrix, mean is the mean value, SD is the standard deviation, control is a normal sample and cancer is a diseased sample.  
         [0055]     In some embodiments, feature selection may be suitable for all classifiers because the categorical characteristics explored and utilized by distinct methods could be quite different  
         [0056]     Alternatively, for a BMM, each parameter may be estimated as the mean of its posterior distribution using, for example, Markov Chain Monte Carlo (MCM) methods. The following will specify the distributions for mutually independent parameters in a BNN model and a GP model, respectively. These two fully Bayesian methods, together with the MLE method for the GP model, will be examined on the National Cancer Institute&#39;s (NCI) ovarian data described in detail below.  
         [0057]     In one embodiment, a mathematically convenient prior of unknown variance is an inverse gamma distribution, the conjugate family for Gaussian weights. In one respect, the conjugate family for Gaussian weights are:
 
 W|σ   0   2   ˜N   u (0, σ 0   2   I   u )  Eq. 12
 
σ 0   2 ˜Inv−Gamma( a   0   , b   0 )  (Eq. 13
 
 M.   .i   |Σ˜N   m (0, Σ)  KEq. 14
 
σ i   2 ˜Inv−Gamma( a   i   , b   i )  Eq. 15
 
 where Inv−Gamma  
           (       θ   ❘   a     ,   b     )     =         b   a       Γ   ⁡     (   a   )         ⁢     θ     -     (     a   +   1     )         ⁢     exp   ⁡     (       -   b     /   θ     )           ,     θ   &gt;     0   ⁢           ⁢   with   ⁢           ⁢   shape   ⁢           ⁢   a     &lt;   0         
 
 and scale b&gt;0. The full conditionals of σ 0   2  and σ i   2  are
 
σ 0   2   |W, M, Σ, D˜ Inv−Gamma(α 0   +u/ 2,  b   0   +∥W∥   2 /2)  Eq. 16
 
σ i   2 |σ 0   2 , . . . , σ i−1   2 , σ i+1   2 , . . . , σ m   2   , W, M, D˜ Inv−Gamma( a   i   +u/ 2,  b   i   +∥M   i ,.∥ 2 /2)  Eq. 17
 
 where a 0 =0.01, b 0 =0.5, a i =0.01, b i =Pr(f i |a) to make the density flatter for larger Pr(f|a). 
 
         [0058]     In a fully Bayesian model of GP, the prior probability of linear part γ 2 =(γ 1   2 , . . . , γ m   2 ) T  and relevance λ 2 =(λ 1   2 , . . . , λ m   2 ) T  are similar to that of M j (Eq. 14), and the prior probability of scale β 2  is similar to that of W i  (Eq. 12). The jitter P and constant cshown in Eq. 6 are set to be 0.01 and 1 respectively.  
       Experiments and Analysis  
       [0059]     For each k-fold cross validation of an embodiment of the MLE-based GP method (k=2, . . . , 10), independent experiments have been repeated 1,000 times on the reduced data by SFS to get the average specificity and sensitivity with a corresponding standard deviation (see  FIG. 4 ). The number of selected features may be controlled by both α and the threshold, for instance on the original binned data, as shown in  FIG. 8 .  
         [0060]     Similarly, MLE-based GP model (GPMLE) generally behaves very stably in leave-one-out cross validations on distinct reduced datasets (see the last row of  FIG. 7 ).  
         [0061]     In one embodiment, a Hybrid Monte Carlo (HMC) method may be used to compute tractable Bayesian integrals by using dynamical methods and Metropolis acceptance rules to propose and accept transition states. In one respect, an NN with one hidden layer of 15 units, where the hidden-to-output weights converge after 80 iterations was analyzed, as shown in  FIGS. 9 and 10 . In  FIG. 9 , the left figure is the HMC simulation of hidden-to-output weights, where u=15 and the z-axis is the number of iterations; in  FIG. 10 , the simulation of the first five input-to-hidden weights of hidden unit  1  are shown. Comparing the GPMLE and SVM, the fully Bayesian methods are both accurate and efficient.  
         [0062]     Although other classifiers may be used, others that were tested were found to be inferior to those by BNN (see  FIG. 11 ). Some sophisticated techniques of data reduction, such as, but not limited to, restriction of coefficient of variation, wavelet analysis, or SVM with soft margin, may also provide good performance of sensitivity and specificity on NCI&#39;s ovarian data. Regardless of the computational complexity, the bagging of MLP is able to yield another reliable Pr(t nεw |x nεw , D), especially for the small training sample, from the viewpoint of ensemble learning. Several samples were found almost always misclassified in cross validations; for example, those explored by MLP bagging in  FIG. 11 .  
         [0063]     Features of the present disclosure may be substantially automated, i.e., computerized. An exemplary computer code for implementing all or parts of this disclosure is shown in  FIG. 12 . The code may be stored in or on any suitable computer readable media, and may be accessed and executed by any suitable computing device, such as a personal desktop or portable computer.  
         [0064]     Although the invention has been disclosed with reference to various particular embodiments, it is understood that equivalents may be employed and substitutions made herein without departing from the scope of the invention as recited in the claims.