/***********************************************************************
Leah Jager
February 2003

This file contains the functions to calculate nonparametric likelihood
confidence bands as described in Owen (1995): 

Owen, A. B.  (1995).  Nonparametric likelihood confidence bands for a 
     distribution function.  Journal of the American Statistical Association 
     90, 516-521.

Functions allow these bands to be calculate in 2 different ways:

(1) By inversion of the Berk-Jones statistic using confidence bands defined
    slightly differently than in Owen and noted in Jager-Wellner
(2) By inversion of the reversed Berk-Jones statistic using confidence bands
    defined in Jager-Wellner
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/**********************************************************************/
double k(double x, double y)
{
  return x*log(x/y) + (1.0-x)*log((1.0-x)/(1.0-y));
}


/**********************************************************************/
double EmpDF(double x, int n, double *data)
     /*  empirical distribution function
	    x is the number at which the function is to be evaluated
	    n is the number of data points
	    data is an array containing the data  */
{
  double tempF_n;
  double F_n;
  int i;

  tempF_n=0.0;
  for (i=1; i <= n; i++) {
    if (data[i] <= x)
      tempF_n += 1;
  }

  F_n = (tempF_n/n);

  return F_n;
}

/*********************************************************************/
double K(double p, double lambda, double x, int n, double data[n])
     /*  function solved to compute H and L for the BJ statistic
	    p is the value we are max/min with respect to
	    lambda is the bound on K
	    x is the number at which we are evaluating the EDF
	    n is the number of data points
	    data is an array containing the data  */
{
  double F_n;
  double k;
  
  F_n = EmpDF(x,n,data);
  k = F_n*log(F_n/p)+(1.0-F_n)*log((1.0-F_n)/(1.0-p))-lambda;
  
  return k;
}

/*********************************************************************/
double K_rev(double p, double lambda, double x, int n, double data[n])
     /*  function solved to compute H and L for the reversed statistic
	    p is the value we are max/min with respect to
	    lambda is the bound on K
	    x is the number at which we are evaluating the EDF
	    n is the number of data points
	    data is an array containing the data  */
{
  double F_n;
  double k;

  F_n = EmpDF(x,n,data);

  k = p*log(p/F_n)+(1.0-p)*log((1.0-p)/(1.0-F_n))-lambda;
  return k;
}

/*********************************************************************/
double F(double x)
     /*  distribution function used in Noe recursion to see with what 
	 probability it is contained in the confidence bands
            x is the number at which the function is evaluated  */
{

  /*x^alpha  alpha=1.0 is uniform*/
  return pow(x,1.0);
  
  /*extreme F1
  return 1.0/(1.0 + log(1.0/x));
  */
  /*extreme F2
  return exp(-((1.0/x) - 1.0));
  */
}

/*********************************************************************/
double PofR_n_BJ(double lambda, double nuisance, double alpha, int n, double *data)
     /*  probability BJ stat is less than lambda using Jager-Wellner a,b
            lambda is lambda
	    nuisance1 is a nuisance parameter
	    alpha is the confidence level
	    n is the number of observations
	    data is the data (should be of length n)  */
{
  int i;
  double *a;
  double *b;
  double answer;
  double x;
  double LL,UL;
  double temp;
  
  a = get_vector(n+1);
  b = get_vector(n+1);
  
  if( !a || !b )
    fprintf(stderr, "Not enough memory to handle n=%d", n);

  b[1]=zbrent(K,lambda,data[0],n,data,1e-16,1.0-1e-16,1e-20);

  for (i=2; i<=n; i++) {

    x = data[i-1];
    
    LL= EmpDF(x,n,data);
    UL= 1.0-1e-16;

    b[i]=zbrent(K,lambda,x,n,data,LL,UL,1e-20);			 
    }
    
  for (i=1; i<=n; i++){
    a[i]=1-b[n-i+1];
   
  }
  
  answer=ecdf_cover(n,a,b,&F)-(1.0-alpha);
 
 /* print out the confidence interval
  printf("a <- c(");
  for(i=1; i<n; i++) {
    printf("%lf,",a[i]);
  }
  printf("%lf)",a[n]);
  printf("b <- c(");
  for(i=1; i<n; i++) {
    printf("%lf,",b[i]);
  }
  printf("%lf)",b[n]);
  */
    /*print out the lambda answer  
      printf("\n\n the answer is %lf", answer);
    */  
  
  free_vector(a, n+1);
  free_vector(b, n+1);
  return answer;
}


/*********************************************************************/
double PofR_n_rev(double lambda, double n_ones, double alpha, int n, double *data)
/*  probability reverse BJ stat is less than lambda using Jager-Wellner a,b
            lambda is lambda
	    n_ones is the number of 1's in b (upper CI)
	    alpha is the confidence level
	    n is the number of observations
	    data is the data (should be of length n)  */
{
  int i,j;
  double *a;
  double *b;
  double answer;
  double x;
  double LL,UL;
  double temp;

  a = get_vector(n+1);
  b = get_vector(n+1);

  if( !a || !b )
    fprintf(stderr, "Not enough memory to handle n=%d", n);

  /* the case for the first order statistic */
  x=data[1];
  b[1]=zbrent(K_rev,lambda,x,n,data,1.0/n,1.0,1e-20);

  int z;   /* the number of ones in the upper confidence bound */  
  z= (int) n_ones;
  j=n-z;

  /* the cases for the other order statistics */
  for (i=2; i<j; i++) {
    
    x = data[i-1];
    LL= EmpDF(x,n,data);
    UL= 1.0;
    b[i]=zbrent(K_rev,lambda,x,n,data,LL,UL,1e-20);
    
}
  x=data[i-1];
  LL=EmpDF(x,n,data);
  UL=1.0;
  b[j]=zbrent(K_rev,lambda,x,n,data,LL,UL,1e-20);

  for(i=(j+1) ;i<=n; i++){
    b[i]=1.0;
  }
  
  for (i=1; i<=n; i++){
    a[i]=1-b[n-i+1];
    
  }
  
  answer=ecdf_cover(n,a,b,&F)-(1.0-alpha);
  /*prints out confidence interval
  printf("a <- c(");
  for(i=1; i<n; i++) {
    printf("%lf,",a[i]);
  }
  printf("%lf)",a[n]);
  printf("b <- c(");
  for(i=1; i<n; i++) {
    printf("%lf,",b[i]);
  }
  printf("%lf)",b[n]);
  */
  /*print out the lambda answer 
    printf("\n\n the answer is %lf", answer);*/

  free_vector(a, n+1);
  free_vector(b, n+1);
  return answer;
  
}
/*********************************************************************/