/*
 * CoxProcess.java
 *
 * Created on 12. Mai 2006, 11:02
 */

package WT;

import java.util.*;


/** Erzeugt ein Realisierungen eines Cox-Prozesses im Intervall [0,T] mit bestimmter Intensitaetsfunktion (sh. Uebungsblatt)
 *
 * @author  Stefanie Eckel
 */
public class CoxProcess {
        
    public static void main(String[] args) {
        
        //Intervallgrenze
        double T = 4;
        
        //Parameter des Intensitaetsprozesses
        double a = 0.5;
        double b = 4;
        double t0 = 2;
        
        //E-Wert
        double e = 0; 
        
        //Leerwahrscheinlichkeit im Intervall [1,2]
        double p = 0;
        
        int iterations = 100;
        for(int i = 0;i<iterations;i++)
        {
            //Erzeugen von auf [0,a] und [0,b] gleichverteilten ZV'en
            double[] U = random(a,b);
            //System.out.println("U1="+U[0]+" U2="+U[1]);
            
            //Erzeugen einer Realisierung des Intensitaetsprozesses und berechnen des Integrals = g(T)
            double gT = IntensityProcess(U, a, b, t0, T);
            //System.out.println("g(t)="+gT);
            
            //g(t0) ausrechnen
            double gt0 = IntensityProcess(U, a, b, t0, t0);
            
            //Erzeugen einer Realisierung eines Poisson-Prozesses und ausgeben der Erneuerungszeitpunkte im Intervall [0,g(T)]
            Vector v = PoissonProcess(1, gT);
            
            //Ruecktransformation der Erneuerungszeitpunkte
            Vector w = transform(U,gt0,T,v);
            
            //Ausgabe einer Realisierung
            /*
            for(int j=0;j<w.size();j++)        
            {
                double s = ((Double) w.get(j)).doubleValue();
                System.out.println(s);
            }
            */
            //Ueberprueft, ob im Intervall [0.1,0.12] ein Erneuerungszeitpunkt liegt
            if(voidProbability(w, 0.1, 0.2)) 
                p++;
            
            //E-Wert
            e = e + w.size();
        }
        
        System.out.println("Erwartungswert: "+e/iterations);
        System.out.println("Leerwahrscheinlichkeit: "+(1-p/iterations));
        
    }
    
    /* Erzeugt zwei ZVen, die auf [0,a] und [0,b] gleichverteilt sind
     *@param a  erste Intervallgrenze
     *@param b  zweite Intervallgrenze
     *
     *@return zwei gleichverteilte ZVen
     */
    public static double[] random(double a, double b)
    {
        //Gleichverteilt auf [0,1]
        double U1 = java.lang.Math.random();
        double U2 = java.lang.Math.random();
        
        //Transformation auf Intervall [0, a] bzw. [0,b]: zur Erinnerung [c,d] => c+(d-c)*U
        U1 = a*U1;
        U2 = b*U2;
           
        return new double[]{U1,U2};
    }
    
    /*Berechnet das Integral ueber eine Realisierung des Intensitaetsprozesses
     @@param U enthaelt die auf [0,a] und [0,b] gleichverteilten ZVen
     *@param a Parameter des Intensitaetsprozesses
     *@param b Parameter des Intensitaetsprozesses
     *@param t0 Parameter des Intensitaetsprozesses
     *@param T Intervallgrenze   
     *
     *@return Wert des Integrals
     */
    public static double IntensityProcess(double[] U, double a, double b, double t0, double T)
    {        
        //U1 bzw. U2
        //Integral analytisch berechnet: U1*t0  bzw. U2*(T-t0)
        double integral = U[0]*t0;
        integral = integral + U[1]*(T-t0);
        
        return integral;
    }
    
    
    /**Erzeugt eine Realisierung des Poisson Prozesses mit Intensitate lambda im Intervall [0,T]
     *
     *@param lambda Intensitaet 
     *@param T Intervallgrenze
     *
     *@return eine Realisierung
     */
    public static Vector PoissonProcess(double lambda, double T)
    {
        double sum = 0;  
        Vector v = new Vector(); 
        
        while(sum<T)
            {
                
                double u = java.lang.Math.random(); //generates a random number between 0 and 1
                double t = - Math.log(u)/lambda; //transforms u into an exponential distributed random variable with intensity lambda
                if(sum+t<T)
                {
                    sum = sum + t;
                    v.add(new Double(sum)); //v.add(double) geht nicht
                    
                    //Ausgabe einer Realisierung
                    //System.out.println(v.size() +" "+sum);
                } 
                else
                    break;
            }      
        return v;
    }
    
    /*Transformiert die Erneuerungszeitpunkte im Intervall [0,g(T)] auf das Intervall [0,T]
     *@param U die zwei gleichverteilten ZVen des Intensitaetsprozesses
     *@param gt0 Integral ueber den Intensitaetsprozesses bis t0
     *@param T Intervallobergrenze
     *@param v Vektor mit den Erneuerungszeitpunkten
     *
     *@return Vektor mit transformierten Erneuerungszeitpunkten
     */
    public static Vector transform(double[] U, double gt0, double T, Vector v)
    {
        Vector w = new Vector();
        
        for(int i=0; i<v.size();i++)
        {
            double s = ((Double) v.get(i)).doubleValue();
            //System.out.println(s);
            double t = 0;
            //g(t)=s
            if(s<gt0) //s < g(t0)
                t = s/U[0];
            else 
                t = T-s/U[1];
            //System.out.println(t);
            w.add(new Double(t));
        }
        return w;
    }
    
    /*Ueberprueft, ob es Erneuerungszeitpunkte im Intervall [a,b] gibt 
     *@param v der Vektor mit Erneuerungszeitpunkten
     *@param a untere Intervallgrenze
     *@param b obere Intervallgrenze
     *
     *@return true falls eine Erneuerungszeitpunkt im Intervall liegt, sonst false
     */
    public static boolean voidProbability(Vector v, double a, double b)
    {
        for(int i=0; i<v.size();i++)
        {    
            double t = ((Double) v.get(i)).doubleValue();
            if(t>=a && t<=b)
            {    
                return true;
            }    
        }
        return false;
    }
}

/*

Ausgabe:
Erwartungswert: 4.5145
Leerwahrscheinlichkeit: 0.7827

*/