Program 6A: Predictions

/*
*/

#ifndef ERROR_LOG_H
#define ERROR_LOG_H

#include <string>

class error_log
{
  public:void check_error (bool should_be_false, const std::string & message);
  void log_error (const std::string & message);
  void clear_error_flag (void);
  void set_error_flag (void);
  bool error_flag (void) const;
    error_log (void);

    protected:bool m_error_flag;
};

#endif

/*
*/

/*
*/

#include "error_log.h"

#include <string>

void
error_log::check_error (bool should_be_false, const std::string & message)
{
  if (should_be_false)
    {
      m_error_flag = true;
      log_error (message);
    }
}

void
error_log::log_error (const std::string & message)
{
  cerr << "Error: " << message << "\n";
}

void
error_log::clear_error_flag (void)
{
  m_error_flag = false;
}

void
error_log::set_error_flag (void)
{
  m_error_flag = true;
}

bool error_log::error_flag (void) const
{
  return m_error_flag;
}

error_log::error_log (void)
{
  clear_error_flag ();
}

/*
*/

#ifndef GAMMA_FUNCTION_H
#define GAMMA_FUNCTION_H

#ifndef SINGLE_VARIABLE_FUNCTION_H
#include "single_variable_function.h"
#endif

class gamma_function:public single_variable_function
{
  public:virtual double at (double x) const;
};

#endif

/*
*/

#include "gamma_function.h"

#ifndef IS_DOUBLE_EQUAL_H
#include "is_double_equal.h"
#endif

#include <math.h>

double
gamma_function::at (double x) const
{
  double
    Result = 0;
  if (is_double_equal (x, 1.0))
    {
      Result = 1;
    }
  else if (is_double_equal (x, 0.5))
    {
      Result = sqrt (M_PI);
    }
  else
    {
      Result = (x - 1.0) * at (x - 1.0);
    }
  return Result;
}

/*
*/

#ifndef IS_DOUBLE_EQUAL_H
#define IS_DOUBLE_EQUAL_H

bool is_double_equal (const double lhs, const double rhs);

void set_is_double_equal_margin (const double new_margin);

#endif

/*
*/

#include "is_double_equal.h"

#include <math.h>

static double double_equal_margin = 0.00001;

bool is_double_equal (const double lhs, const double rhs)
{
  if (fabs (lhs - rhs) < double_equal_margin)
    {
      return true;
    }
  else
    {
      return false;
    }
}

void
set_double_equal_margin (const double new_margin)
{
  double_equal_margin = new_margin;
}

/*
*/

#ifndef PAIRED_NUMBER_LIST_PREDICTOR_H
#define PAIRED_NUMBER_LIST_PREDICTOR_H

#ifndef PAIRED_NUMBER_LIST_H
#include "paired_number_list.h"
#endif
#ifndef T_DISTRIBUTIONL_H
#include "t_distribution.h"
#endif

class paired_number_list_predictor:public paired_number_list
{
  public:double variance (void) const;
  double standard_deviation (void) const;
  double projected_y (double x) const;
  double prediction_range (double x, double range) const;
  double lower_prediction_interval (double x, double range) const;
  double upper_prediction_interval (double x, double range) const;
  double t (double range) const;

    protected:t_distribution m_t_distribution;
  double prediction_range_base (void) const;
};

#endif

/*
*/

#include "paired_number_list_predictor.h"

#ifndef CONTRACT_H
#include "contract.h"
#endif

#include <math.h>

double
paired_number_list_predictor::variance (void) const
{
  REQUIRE (entry_count () > 2);
  double
    Result = 0;
  list < double >::const_iterator x_iter;
  list < double >::const_iterator y_iter;

  for (x_iter = m_xs.begin (), y_iter = m_ys.begin ();
       (x_iter != m_xs.end ()) && (y_iter != m_ys.end ()); ++x_iter, ++y_iter)
    {
      Result += pow (*y_iter - beta_0 () - beta_1 () * (*x_iter), 2);
    }
  Result *= 1.0 / (entry_count () - 2.0);
  return Result;
}

double
paired_number_list_predictor::standard_deviation (void) const
{
  return sqrt (variance ());
}

double
paired_number_list_predictor::projected_y (double x) const
{
  return beta_0 () + beta_1 () * x;
}

double
paired_number_list_predictor::t (double range) const
{
  const_cast <
    paired_number_list_predictor *
    >(this)->m_t_distribution.set_n (entry_count () - 2);
  return m_t_distribution.at (range);
}

double
paired_number_list_predictor::prediction_range (double x, double range) const
{
  REQUIRE (entry_count () > 0);
  const double
    a_t = t (range);
  const double
    dev = standard_deviation ();
  const double
    x_m = x_mean ();
  double
    ecount_inv = 1 / entry_count ();
  return t (range) * standard_deviation ()
    * sqrt (1.0 + 1.0 / static_cast < double >(entry_count ())
	    + pow (x - x_mean (), 2) / prediction_range_base ());
}

double
paired_number_list_predictor::lower_prediction_interval (double x,
							 double range) const
{
  return projected_y (x) - prediction_range (x, range);
}

double
paired_number_list_predictor::upper_prediction_interval (double x,
							 double range) const
{
  return projected_y (x) + prediction_range (x, range);
}

double
paired_number_list_predictor::prediction_range_base (void) const
{
  double
    Result = 0;
  for (std::list < double >::const_iterator x_iter = m_xs.begin ();
       x_iter != m_xs.end (); ++x_iter)
    {
      Result += pow ((*x_iter) - x_mean (), 2);
    }
  return Result;
}

/*
*/

#ifndef T_DISTRIBUTION_H
#define T_DISTRIBUTION_H

#ifndef SINGLE_VARIABLE_FUNCTION_H
#include "single_variable_function.h"
#endif

class t_distribution:public single_variable_function
{
  public:virtual double at (double x) const;
  void set_n (double new_n);
    t_distribution (void);

    protected:double n;
  double next_guess (double arg, double last_result, double target) const;
};

#endif

#include "t_distribution.h"
#ifndef T_INTEGRAL_H
#include "t_integral.h"
#endif

void
t_distribution::set_n (double new_n)
{
  n = new_n;
}

double
t_distribution::at (double x) const
{
  //x is given in a range; we need to convert this to a two-sided
  //distribution so we get what we're expecting
  x = 0.5 + x / 2;
  t_integral
    integral;
  integral.set_n (n);
  double
    last_error = 0;
  const double
    error_margin = 0.0000000001;
  double
    integral_arg = 0;
  double
    this_result = 0.5;
  double
    last_result = 0;
  bool
    has_tried_once = false;
  while (!has_tried_once || (last_error > error_margin))
    {
      last_result = this_result;
      this_result = integral.at (integral_arg);
      last_error = fabs (this_result - x);
      //store the argument in a temp (don't assign to last_arg,
      //because we need both values for the next guess!
      integral_arg = next_guess (integral_arg, last_result, x);
      has_tried_once = true;
    }
  return integral_arg;
}

t_distribution::t_distribution (void)
{
  set_n (0);
}

double
t_distribution::next_guess (double arg, double last_result, double target) const
{
  double
    Result = 0;
  Result = arg + (target - last_result);
  return Result;
}

#ifndef T_DISTRIBUTION_BASE_H
#define T_DISTRIBUTION_BASE_H

#ifndef SINGLE_VARIABLE_FUNCTION_H
#include "single_variable_function.h"
#endif

class t_distribution_base:public single_variable_function
{
  public:virtual double at (double x) const;
  void set_n (double new_n);
    t_distribution_base (void);
    protected:double n;
};

#endif

/*
*/

#include "t_distribution_base.h"
#ifndef CONTRACT_H
#include "contract.h"
#endif
#include <math.h>


double
t_distribution_base::at (double x) const
{
  REQUIRE (n > 0);
  return pow (1 + x * x / n, -(n + 1) / 2);
}

void
t_distribution_base::set_n (double new_n)
{
  n = new_n;
}

t_distribution_base::t_distribution_base (void):
n (0)
{
}


/*
*/

#ifndef T_INTEGRAL_H
#define T_INTEGRAL_H

#ifndef SINGLE_VARIABLE_FUNCTION_H
#include "single_variable_function.h"
#endif
#ifndef GAMMA_FUNCTION_H
#include "gamma_function.h"
#endif
#ifndef SIMPSON_INTEGRATOR_H
#include "simpson_integrator.h"
#endif
#ifndef T_DISTRIBUTION_BASE_H
#include "t_distribution_base.h"
#endif

class t_integral:public single_variable_function
{
  public:virtual double at (double x) const;
  void set_n (double new_n);
    t_integral (void);

    protected:simpson_integrator homer;
  gamma_function gamma;
  t_distribution_base base;
  double n;
  double multiplier (void) const;
};

#endif

/*
*/

#include "t_integral.h"
#ifndef CONTRACT_H
#include "contract.h"
#endif

#include <math.h>


double
t_integral::multiplier (void) const
{
  return gamma.at ((n + 1) / 2) / (sqrt (n * M_PI) * gamma.at (n / 2));
}

void
t_integral::set_n (double new_n)
{
  n = new_n;
  base.set_n (new_n);
}

t_integral::t_integral (void)
{
  set_n (0);
}

double
t_integral::at (double x) const
{
  REQUIRE (n > 0);
  double
    Result = 0.5;
  if (x < 0)
    {
      Result = 0.5 - multiplier () * homer.integral (base, 0, -x);
    }
  else if (x > 0)
    {
      Result = 0.5 + multiplier () * homer.integral (base, 0, x);
    }
  return Result;
}

/*
*/

/*
*/

#include <fstream>
#include <iostream>
#include "string.h"

#ifndef PREDICTOR_PARSER_H
#include "predictor_parser.h"
#endif

istream *
input_stream_from_args (int arg_count, const char **arg_vector)
{
  istream *Result = NULL;
  if (arg_count == 1)
    {
      Result = &cin;
    }
  else
    {
      const char *help_text =
	"PSP exercise 6A: Calculate a prediction and interval given historical data.\nUsage:\n\tpsp_5a\n\n";
      cout << help_text;
    }
  return Result;
}

int
main (int arg_count, const char **arg_vector)
{
  //get the input stream, or print the help text as appropriate
  istream *input_stream = input_stream_from_args (arg_count, arg_vector);
  if (input_stream != NULL)
    {
      predictor_parser parser;
      parser.set_input_stream (input_stream);
      parser.parse_until_eof ();
    }
}

/*
*/

class ERROR_LOG
--very rudimentary error reporting and error flag tracking.

creation {ANY} 
   make

feature {ANY} 
   
   check_for_error(should_be_false: BOOLEAN; message: STRING) is 
      do  
         if should_be_false then 
            log_error(message);
         end; 
      end -- check_for_error
   
   log_error(message: STRING) is 
      do  
         std_error.put_string(message + "%N");
         set_error_flag;
      end -- log_error
   
   clear_error_flag is 
      do  
         error_flag := false;
      end -- clear_error_flag
   
   set_error_flag is 
      do  
         error_flag := true;
      end -- set_error_flag
   
   error_flag: BOOLEAN;
   
   make is 
      do  
         clear_error_flag;
      end -- make

end -- class ERROR_LOG

class GAMMA_FUNCTION
--gamma function, one of the bases of the t-distribution

inherit 
   SINGLE_VARIABLE_FUNCTION
      redefine at
      end; 
   
feature {ANY} 
   
   error_margin: DOUBLE is 0.000001;
   
   acceptably_equal(lhs, rhs: DOUBLE): BOOLEAN is 
      --is lhs within error_margin of rhs?
      do  
         Result := false;
         if lhs.in_range(rhs - error_margin,rhs + error_margin) then 
            Result := true;
         end; 
      end -- acceptably_equal
   
   at(x: DOUBLE): DOUBLE is 
      --gamma function
      do  
         if acceptably_equal(x,1.0) then 
            Result := 1;
         elseif acceptably_equal(x,0.5) then 
            Result := Pi.sqrt;
         else 
            Result := (x - 1.0) * at(x - 1.0);
         end; 
      end -- at

end -- class GAMMA_FUNCTION

class MAIN

creation {ANY} 
   make

feature {ANY} 
   
   make is 
      local 
         parser: PREDICTOR_PARSER;
         gamma: GAMMA_FUNCTION;
      do  
         !!parser.make;
         parser.set_input(io);
         parser.parse_until_eof;
      end -- make

end -- MAIN

class PAIRED_NUMBER_LIST_PREDICTOR
--reads a set of paired numbers, does linear regression, predicts results

inherit 
   PAIRED_NUMBER_LIST
      redefine make
      end; 
   
creation {ANY} 
   make

feature {ANY} 
   
   variance: DOUBLE is 
      local 
         i: INTEGER;
      do  
         Result := 0;
         from 
            i := xs.lower;
         until 
            not (xs.valid_index(i) and ys.valid_index(i))
         loop 
            Result := Result + (ys.item(i) - beta_0 - beta_1 * xs.item(i)) ^ 2;
            i := i + 1;
         end; 
         Result := Result / (entry_count - 2);
      end -- variance
   
   standard_deviation: DOUBLE is 
      do  
         Result := variance.sqrt;
      end -- standard_deviation
   
   projected_y(x: DOUBLE): DOUBLE is 
      --projected value of given x, using linear regression 
      --parameters from xs and ys
      do  
         Result := beta_0 + beta_1 * x;
      end -- projected_y
   
   prediction_range_base: DOUBLE is 
      --base of the prediction range, used in prediction_range
      local 
         i: INTEGER;
      do  
         Result := 0;
         from 
            i := xs.lower;
         until 
            not (xs.valid_index(i) and ys.valid_index(i))
         loop 
            Result := Result + (xs.item(i) - xs.mean) ^ 2;
            i := i + 1;
         end; 
      end -- prediction_range_base
   
   prediction_range(x, range: DOUBLE): DOUBLE is 
      --prediction range, based on given estimate and % range
      require 
         entry_count > 0; 
      do  
         Result := (1.0 + (1.0 / entry_count.to_double) + (((x - xs.mean) ^ 2) / prediction_range_base)).sqrt;
         Result := t(range) * standard_deviation * Result;
      end -- prediction_range
   
   lower_prediction_interval(x, range: DOUBLE): DOUBLE is 
      --LPI, from [Humphrey95]
      do  
         Result := projected_y(x) - prediction_range(x,range);
      end -- lower_prediction_interval
   
   upper_prediction_interval(x, range: DOUBLE): DOUBLE is 
      --UPI, from [Humphrey95]
      do  
         Result := projected_y(x) + prediction_range(x,range);
      end -- upper_prediction_interval
   
   t_distribution: T_DISTRIBUTION;
   
   make is 
      do  
         Precursor;
         !!t_distribution.make;
      end -- make
   
   t(range: DOUBLE): DOUBLE is 
      --gets the size of the t-distribution at the given alpha range
      do  
         t_distribution.set_n(entry_count - 2);
         Result := t_distribution.at(range);
      end -- t

end -- class PAIRED_NUMBER_LIST_PREDICTOR

class PREDICTOR_PARSER
--reads a list of number pairs, and performs linear regression analysis

inherit 
   SIMPLE_INPUT_PARSER
      redefine parse_last_line, transformed_line
      end; 
   
creation {ANY} 
   make

feature {ANY} 
   
   inline_comment_begin: STRING is "--";
   
   string_stripped_of_comment(string: STRING): STRING is 
      --strip the string of any comment
      local 
         comment_index: INTEGER;
      do  
         if string.has_string(inline_comment_begin) then 
            comment_index := string.index_of_string(inline_comment_begin);
            if comment_index = 1 then 
               Result := "";
            else 
               Result := string.substring(1,comment_index - 1);
            end; 
         else 
            Result := string;
         end; 
      end -- string_stripped_of_comment
   
   string_stripped_of_whitespace(string: STRING): STRING is 
      --strip string of whitespace
      do  
         Result := string;
         Result.left_adjust;
         Result.right_adjust;
      end -- string_stripped_of_whitespace
   
   transformed_line(string: STRING): STRING is 
      --strip comments and whitespace from parseable line      
      do  
         Result := string_stripped_of_whitespace(string_stripped_of_comment(string));
      end -- transformed_line
   
   number_list: PAIRED_NUMBER_LIST_PREDICTOR;

feature {ANY} --parsing

   found_end_of_historical_data: BOOLEAN;
   
   reset is 
      --resets the parser and makes it ready to go again
      do  
         found_end_of_historical_data := false;
         number_list.reset;
      end -- reset
   
   make is 
      do  
         !!number_list.make;
         reset;
      end -- make
   
   parse_last_line_as_historical_data is 
      --interpret last_line as a pair of comma-separated values
      local 
         error_log: ERROR_LOG;
         comma_index: INTEGER;
         x_string: STRING;
         y_string: STRING;
         new_x: DOUBLE;
         new_y: DOUBLE;
      do  
         !!error_log.make;
         comma_index := last_line.index_of(',');
         error_log.check_for_error(comma_index = last_line.count + 1,"No comma:" + last_line);
         x_string := last_line.substring(1,comma_index - 1);
         y_string := last_line.substring(comma_index + 1,last_line.count);
         error_log.check_for_error(not (x_string.is_double or x_string.is_integer),"invalid X:" + last_line);
         error_log.check_for_error(not (y_string.is_double or y_string.is_integer),"invalid Y:" + last_line);
         if not error_log.error_flag then 
            new_x := double_from_string(x_string);
            new_y := double_from_string(y_string);
            number_list.add_entry(new_x,new_y);
            std_output.put_string("added: ");
            std_output.put_double(new_x);
            std_output.put_string(", ");
            std_output.put_double(new_y);
            std_output.put_new_line;
         end; 
      end -- parse_last_line_as_historical_data
   
   double_from_string(string: STRING): DOUBLE is 
      require 
         string.is_double or string.is_integer; 
      do  
         if string.is_double then 
            Result := string.to_double;
         elseif string.is_integer then 
            Result := string.to_integer.to_double;
         end; 
      end -- double_from_string
   
   historical_data_terminator: STRING is "stop";
   
   parse_last_line_as_end_of_historical_data is 
      --interpret last line as the end of historical data
      require 
         last_line.compare(historical_data_terminator) = 0; 
      do  
         found_end_of_historical_data := true;
         std_output.put_string("Historical data read.%NBeta-0: ");
         std_output.put_double(number_list.beta_0);
         std_output.put_string("%NBeta-1: ");
         std_output.put_double(number_list.beta_1);
         std_output.put_string("%NStandard Deviation: ");
         std_output.put_double(number_list.standard_deviation);
         std_output.put_string("%N%N");
      end -- parse_last_line_as_end_of_historical_data
   
   parse_last_line_as_prediction is 
      --interpret last line as a single x, for a predictive y
      local 
         error_log: ERROR_LOG;
         x: DOUBLE;
      do  
         !!error_log.make;
         error_log.check_for_error(not (last_line.is_double or last_line.is_integer),"Not a double : " + last_line);
         if not error_log.error_flag then 
            x := double_from_string(last_line);
            std_output.put_string("Estimate at x=");
            std_output.put_double(x);
            std_output.put_string("%NProjected y: ");
            std_output.put_double(number_list.projected_y(x));
            std_output.put_string("%Nt (70 percent): ");
            std_output.put_double(number_list.t(0.7));
            std_output.put_string("%Nt (90 percent): ");
            std_output.put_double(number_list.t(0.9));
            std_output.put_string("%NRange (70 percent): ");
            std_output.put_double(number_list.prediction_range(x,0.7));
            std_output.put_string("; UPI: ");
            std_output.put_double(number_list.upper_prediction_interval(x,0.7));
            std_output.put_string("; LPI: ");
            std_output.put_double(number_list.lower_prediction_interval(x,0.7));
            std_output.put_string("%NRange (90 percent): ");
            std_output.put_double(number_list.prediction_range(x,0.9));
            std_output.put_string("; UPI: ");
            std_output.put_double(number_list.upper_prediction_interval(x,0.9));
            std_output.put_string("; LPI: ");
            std_output.put_double(number_list.lower_prediction_interval(x,0.9));
            std_output.put_new_line;
         end; 
      end -- parse_last_line_as_prediction
   
   parse_last_line is 
      --parse the last line according to state
      do  
         if not last_line.empty then 
            if last_line.compare(historical_data_terminator) = 0 then 
               parse_last_line_as_end_of_historical_data;
            else 
               if found_end_of_historical_data then 
                  parse_last_line_as_prediction;
               else 
                  parse_last_line_as_historical_data;
               end; 
            end; 
         end; 
      end -- parse_last_line

end -- class PREDICTOR_PARSER

class T_DISTRIBUTION
--the t-distribution, used to find prediction ranges, etc.

inherit 
   SINGLE_VARIABLE_FUNCTION
      redefine at
      end; 
   
creation {ANY} 
   make

feature {ANY} 
   
   n: INTEGER;
   
   set_n(new_n: INTEGER) is 
      do  
         n := new_n;
      end -- set_n
   
   make is 
      do  
         set_n(0);
      end -- make
   
   next_guess(arg, last_result, target: DOUBLE): DOUBLE is 
      do  
         Result := arg + (target - last_result);
      end -- next_guess
   
   error_margin: DOUBLE is 0.0000000001;
   
   at(x: DOUBLE): DOUBLE is 
      local 
         t_integral: T_INTEGRAL;
         last_error: DOUBLE;
         this_result: DOUBLE;
         last_result: DOUBLE;
         has_tried_once: BOOLEAN;
         target: DOUBLE;
      do  
         --convert the range to a two-sided distribution so we get 
         --what we're expecting
         target := 0.5 + x / 2;
         !!t_integral.make;
         t_integral.set_n(n);
         last_error := 0;
         from 
            Result := 0;
            has_tried_once := false;
         until 
            has_tried_once and error_margin > last_error
         loop 
            last_result := this_result;
            this_result := t_integral.at(Result);
            last_error := (this_result - target).abs;
            Result := next_guess(Result,last_result,target);
            has_tried_once := true;
         end; 
      end -- at

end -- class T_DISTRIBUTION

class T_DISTRIBUTION_BASE
--simple function as the base of the t-distribution integral

inherit 
   SINGLE_VARIABLE_FUNCTION
      redefine at
      end; 
   
creation {ANY} 
   make

feature {ANY} 
   
   n: INTEGER;
      --degrees of freedom
   
   set_n(new_n: INTEGER) is 
      do  
         n := new_n;
      end -- set_n
   
   make is 
      do  
         set_n(0);
      end -- make
   
   at(x: DOUBLE): DOUBLE is 
      require 
         n > 0; 
      do  
         Result := ((1.0 + x * x / n.to_double) ^ - (n + 1)).sqrt;
      end -- at

end -- class T_DISTRIBUTION_BASE

class T_INTEGRAL
--integral of the t-distribution up to a given point; NOT what you 
--use to make predictions; use T_DISTRIBUTION instead

inherit 
   SINGLE_VARIABLE_FUNCTION
      redefine at
      end; 
   
creation {ANY} 
   make

feature {ANY} 
   
   n: INTEGER;
   
   set_n(new_n: INTEGER) is 
      do  
         n := new_n;
         t_distribution_base.set_n(new_n);
      end -- set_n
   
   make is 
      do  
         !!gamma;
         !!homer.make;
         !!t_distribution_base.make;
         set_n(0);
      end -- make
   
   multiplier: DOUBLE is 
      do  
         Result := gamma.at((n + 1) / 2) / ((n * Pi).sqrt * gamma.at(n / 2));
      end -- multiplier
   
   gamma: GAMMA_FUNCTION;
   
   homer: SIMPSON_INTEGRATOR;
   
   t_distribution_base: T_DISTRIBUTION_BASE;
   
   at(x: DOUBLE): DOUBLE is 
      do  
         Result := 0.5;
         if x < 0 then 
            Result := 0.5 - multiplier * homer.integral(t_distribution_base,0,- x);
         else 
            Result := 0.5 + multiplier * homer.integral(t_distribution_base,0,x);
         end; 
      end -- at

end -- class T_INTEGRAL