mev-4.0.01/source/org/tigr/microarray/mev/motif/BetaPrior.java

Code
Comments
Other
Rev Date Author Line
2 26 Feb 07 jari 1 /*
2 26 Feb 07 jari 2 Copyright @ 1999-2005, The Institute for Genomic Research (TIGR).
2 26 Feb 07 jari 3 All rights reserved.
2 26 Feb 07 jari 4 */
2 26 Feb 07 jari 5 /*
2 26 Feb 07 jari 6  * $RCSfile: BetaPrior.java,v $
2 26 Feb 07 jari 7  * $Revision: 1.3 $
2 26 Feb 07 jari 8  * $Date: 2005/03/10 15:40:35 $
2 26 Feb 07 jari 9  * $Author: braistedj $
2 26 Feb 07 jari 10  * $State: Exp $
2 26 Feb 07 jari 11  */
2 26 Feb 07 jari 12 package org.tigr.microarray.mev.motif;
2 26 Feb 07 jari 13
2 26 Feb 07 jari 14 /*************************************************************************
2 26 Feb 07 jari 15  *
2 26 Feb 07 jari 16  * beta prior distribution for Bayesean statistics:
2 26 Feb 07 jari 17  *
2 26 Feb 07 jari 18  * mean = m = a/(a+b)
2 26 Feb 07 jari 19  * var = S^2 = a*b/((a+b+1)(a+b)^2)
2 26 Feb 07 jari 20  *
2 26 Feb 07 jari 21  * a = -m - m^2*(m-1)/s^2
2 26 Feb 07 jari 22  * b = -1 +m + m*(m-1)^2/s^2
2 26 Feb 07 jari 23  *
2 26 Feb 07 jari 24  *************************************************************************/
2 26 Feb 07 jari 25
2 26 Feb 07 jari 26 public class BetaPrior {
2 26 Feb 07 jari 27     
2 26 Feb 07 jari 28     double A;        /* total pseudo trials */
2 26 Feb 07 jari 29     double N;        /* total real trials */
2 26 Feb 07 jari 30     double T;        /* Total trials */
2 26 Feb 07 jari 31     double alpha;    /* #successful pseudo trials */
2 26 Feb 07 jari 32     double beta;     /* #failed pseudo trials */
2 26 Feb 07 jari 33     long success;    /* #successful real trials */
2 26 Feb 07 jari 34     long expect;     /* expected number of success */
2 26 Feb 07 jari 35     double p;        /* posterior probability of sucess */
2 26 Feb 07 jari 36     double weight;   /* A/N = w/(1-w) - fractional weight */
2 26 Feb 07 jari 37     boolean calc;
2 26 Feb 07 jari 38     
2 26 Feb 07 jari 39     /*******************************************************************
2 26 Feb 07 jari 40      * N = total sites;   - Pseudo and Real sites.
2 26 Feb 07 jari 41      * tot_sites  = number of sites;
2 26 Feb 07 jari 42      * alpha = number site pseudo counts;
2 26 Feb 07 jari 43      * beta = number site pseudo counts;
2 26 Feb 07 jari 44      *******************************************************************/
2 26 Feb 07 jari 45     public BetaPrior(long Expect, double Weight, double n) {
2 26 Feb 07 jari 46   if (N <=0) System.out.println("total sites must be > 0");
2 26 Feb 07 jari 47   expect = Expect;
2 26 Feb 07 jari 48   weight = Weight;
2 26 Feb 07 jari 49   SetBPriorN(N);
2 26 Feb 07 jari 50   success = 0;
2 26 Feb 07 jari 51   p = (alpha + (double)expect)/T;
2 26 Feb 07 jari 52   calc = false;
2 26 Feb 07 jari 53     }
2 26 Feb 07 jari 54     
2 26 Feb 07 jari 55     
2 26 Feb 07 jari 56     public void ClearBPrior() {
2 26 Feb 07 jari 57   success = 0;
2 26 Feb 07 jari 58   p = (alpha + (double)expect)/T;
2 26 Feb 07 jari 59   calc = false;
2 26 Feb 07 jari 60     }
2 26 Feb 07 jari 61     
2 26 Feb 07 jari 62     public void SetBPriorS(long Success) {
2 26 Feb 07 jari 63   success = Success;
2 26 Feb 07 jari 64   calc = true;
2 26 Feb 07 jari 65     }
2 26 Feb 07 jari 66     
2 26 Feb 07 jari 67     public void SetBPriorN(double n) {
2 26 Feb 07 jari 68   double  ratio;
2 26 Feb 07 jari 69   ratio = (weight/(1.0 - weight));
2 26 Feb 07 jari 70   this.N = n;
2 26 Feb 07 jari 71   A =n*ratio;     /* A = N*(w/(1-w)) */
2 26 Feb 07 jari 72   alpha = (double) expect*ratio;
2 26 Feb 07 jari 73   beta= A - alpha;
2 26 Feb 07 jari 74   T=(A + N);
2 26 Feb 07 jari 75   calc = true;
2 26 Feb 07 jari 76     }
2 26 Feb 07 jari 77     
2 26 Feb 07 jari 78     public double PostProbBPrior() {
2 26 Feb 07 jari 79   if (calc) {
2 26 Feb 07 jari 80       calc=false;
2 26 Feb 07 jari 81       p=(alpha+(double)success)/(T-1);
2 26 Feb 07 jari 82   }
2 26 Feb 07 jari 83   return p;
2 26 Feb 07 jari 84     }
2 26 Feb 07 jari 85     
2 26 Feb 07 jari 86     /** Return the ratio of new to old. **/
2 26 Feb 07 jari 87     public double RatioBPrior(BetaPrior B1, BetaPrior B2) {
2 26 Feb 07 jari 88   double v1,v2;
2 26 Feb 07 jari 89   v1 = (double) B1.success;
2 26 Feb 07 jari 90   v2 = (double) B2.success;
2 26 Feb 07 jari 91   v1 = LnBeta(v1+B1.alpha,B1.N - v1 + B1.beta)- LnBeta(B1.alpha,B1.beta);
2 26 Feb 07 jari 92   v2 = LnBeta(v2+B2.alpha,B2.N - v2 + B2.beta)- LnBeta(B2.alpha,B2.beta);
2 26 Feb 07 jari 93   return Math.exp(v1 - v2);
2 26 Feb 07 jari 94     }
2 26 Feb 07 jari 95     
2 26 Feb 07 jari 96     public double LnBeta(double a, double b) {
2 26 Feb 07 jari 97   return(lgamma((float)(a))+lgamma((float)(b))-lgamma((float)(a+b)));
2 26 Feb 07 jari 98     }
2 26 Feb 07 jari 99     
2 26 Feb 07 jari 100     public float lgamma(float xx) {
2 26 Feb 07 jari 101   double x,y,tmp,ser;
2 26 Feb 07 jari 102   double[] cof=new double[6];
2 26 Feb 07 jari 103   cof[0]=76.18009172947146;
2 26 Feb 07 jari 104   cof[1]=-86.50532032941677;
2 26 Feb 07 jari 105   cof[2]=24.01409824083091;
2 26 Feb 07 jari 106   cof[3]=-1.231739572450155;
2 26 Feb 07 jari 107   cof[4]=0.1208650973866179e-2;
2 26 Feb 07 jari 108   cof[5]=-0.5395239384953e-5;
2 26 Feb 07 jari 109   int j;
2 26 Feb 07 jari 110   y=x=xx;
2 26 Feb 07 jari 111   tmp=x+5.5;
2 26 Feb 07 jari 112   tmp -= (x+0.5)*Math.log(tmp);
2 26 Feb 07 jari 113   ser=1.000000000190015;
2 26 Feb 07 jari 114   for (j=0;j<=5;j++) ser += cof[j]/++y;
2 26 Feb 07 jari 115   return(float)(-tmp+Math.log(2.5066282746310005*ser/x));
2 26 Feb 07 jari 116     }
2 26 Feb 07 jari 117 }

Generated on Thu Jul 20 09:30:03 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)