/*
 * Distributions.java
 *
 * Created on March 8, 2006, 11:45 AM
 *
 * To change this template, choose Tools | Options and locate the template under
 * the Source Creation and Management node. Right-click the template and choose
 * Open. You can then make changes to the template in the Source Editor.
 */

package probabilitycalculator;
import javax.swing.JOptionPane;

/**
 *
 * @author Rinaman
 */
public abstract class Distributions {
    
    /** Creates a new instance of Distributions */
    public Distributions() {
    }
    public double rtflsp(double x, double x1, double x2, double xacc, double parmone, double parmtwo) {
        int MAXIT = 1000;
        int j;
        double fl, fh, xl, xh, swap, dx, del, f, rtf;
        fl = x - df(x1, parmone, parmtwo);
        fh = x - df(x2, parmone, parmtwo);
        if(fl*fh > 0.0) {
            JOptionPane.showMessageDialog(null,"Error: Root not bracketed in rtflsp");
            return -2.0e12;
        }
        if(fl < 0.0) {
            xl = x1;
            xh = x2;
        }
        else {
            xl = x2;
            xh = x1;
            swap = fl;
            fl = fh;
            fh = swap;
        }
        dx = xh - xl;
        for(j = 1;j <= MAXIT;j++) {
            rtf = xl + dx*fl/(fl-fh);
            f = x - df(rtf, parmone, parmtwo);
            if(f < 0.0) {
                del = xl - rtf;
                xl = rtf;
                fl = f;
            }
            else {
                del = xh - rtf;
                xh = rtf;
                fh = f;
            }
            dx = xh - xl;
            if(Math.abs(del) < xacc || f == 0.0)return rtf;
        }
        JOptionPane.showMessageDialog(null,"Error: Maximum iterations exceeded in rtflsp.");
        return -1.0;
    }
    
    public abstract double df(double x, double parmone, double parmtwo);
    
    public double erfc(double x) {
        return x < 0.0 ? gammp(0.5,x*x) : gammq(0.5,x*x);
    }
    
    public double gammp(double a, double x) {
        if(x < (a + 1.0))return gser(a,x);
        else return 1.0 - gcf(a,x);
    }
    
    public double gammq(double a, double x) {
        if(x < (a + 1.0))return 1.0 - gser(a,x);
        else return gcf(a,x);
    }
    
    public double gser(double a, double x) {
        double EPS = 3.0e-7;
        double sum, del, ap, gln;
        int n;
        gln = gammln(a);
        ap = a;
        del = sum = 1.0/a;
        for(n = 1;n <= 1000;n++) {
            ap += 1.0;
            del *= x/ap;
            sum += del;
            if(Math.abs(del) < Math.abs(sum)*EPS)
                return sum*Math.exp(-x+a*Math.log(x)-gln);
        }
        JOptionPane.showMessageDialog(null,"Error: Value out of range in gser");
        return -2.0;
    }
    
    public double gcf(double a, double x) {
        double EPS = 3.0e-7;
        int n;
        double gln, gold = 0.0, g, fac = 1.0, b1 = 1.0;
        double b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
        gln = gammln(a);
        a1 = x;
        for(n = 1;n <= 1000;n++) {
            an = (double)n;
            ana = an - a;
            a0 = (a1+a0*ana)*fac;
            b0 = (b1+b0*ana)*fac;
            anf = an*fac;
            a1 = x*a0+anf*a1;
            b1 = x*b0+anf*b1;
            if(a1 != 0.0) {
                fac = 1.0/a1;
                g = b1*fac;
                if(Math.abs((g-gold)/g) < EPS)
                    return Math.exp(-x+a*Math.log(x)-gln)*g;
                gold = g;
            }
        }
        return -2.0;
    }
    
    public double gammln(double xx) {
        double x, tmp, ser;
        double cof[] = {76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5};
        int j;
        x = xx - 1.0;
        tmp = x + 5.5;
        tmp -= (x+0.5)*Math.log(tmp);
        ser = 1.0;
        for(j = 0;j < 6;j++) {
            x += 1.0;
            ser += cof[j]/x;
        }
        return -tmp + Math.log(2.50662827465*ser);
    }
    
    public double betai(double a, double b, double x) {
        double bt;
        if(x == 0.0 || x == 1.0)bt = 0.0;
        else bt = Math.exp(gammln(a+b)-gammln(a)-gammln(b)+a*Math.log(x)+b*Math.log(1.0-x));
        if(x < (a+1.0)/(a+b+2.0))return bt*betacf(a,b,x)/a;
        else return 1.0 - bt*betacf(b,a,1.0-x)/b;
    }
    
    public double betacf(double a, double b, double x) {
        double EPS = 3.0e-7;
        double qap, qam, qab, em, tem, d;
        double bz, bm = 1.0, bp, bpp;
        double az = 1.0, am = 1.0, ap, app, aold;
        int m;
        qab = a + b;
        qap = a + 1.0;
        qam = a - 1.0;
        bz = 1.0 - qab*x/qap;
        for(m = 1; m < 1000; m++) {
            em = (double)m;
            tem = em + em;
            d = em*(b-em)*x/((qam+tem)*(a+tem));
            ap = az + d*am;
            bp = bz + d*bm;
            d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
            app = ap + d*az;
            bpp = bp + d*bz;
            aold = az;
            am = ap/bpp;
            bm = bp/bpp;
            az = app/bpp;
            bz = 1.0;
            if(Math.abs(az-aold) < (EPS*Math.abs(az)))return az;
        }
        JOptionPane.showMessageDialog(null,"Error: a or b too big of too many iterations in Betacf");
        return 0.0;
    }
}