File Coverage

blib/lib/Statistics/KernelEstimation.pm

Criterion	Covered	Total	%
statement	198	266	74.4
branch	51	104	49.0
condition	17	30	56.6
subroutine	33	38	86.8
pod	17	17	100.0
total	316	455	69.4

package Statistics::KernelEstimation;

use 5.008008;
use strict;
use warnings;

use Carp;

our $VERSION = '0.05';

# =================================================================
# TO DOs
#
# - More unit tests
#   - bandidth from data
#   - math function
#   - optimization
#
# - opt broken for epanechnikov
# - max number of integration steps in stiffness integral

# =================================================================
# Ctors

sub new {
  return new_gauss( @_ );
}

sub new_gauss {
  my $self = _new( @_ );

  $self->{pdf} = \&_gauss_pdf;
  $self->{cdf} = \&_gauss_cdf;

  $self->{curvature} = \&_gauss_curvature;
  $self->{extension} = 3;

  return $self;
}

sub new_box {
  my $self = _new( @_ );

  $self->{pdf} = \&_box_pdf;
  $self->{cdf} = \&_box_cdf;

  $self->{curvature} = \&_box_curvature;
  $self->{extension} = 0;

  $self->{optimizable} = 0;

  return $self;

}

sub new_epanechnikov {
  my $self = _new( @_ );

  $self->{pdf} = \&_epanechnikov_pdf;
  $self->{cdf} = \&_epanechnikov_cdf;

  $self->{curvature} = \&_epanechnikov_curvature;
  $self->{extension} = 0;

  $self->{optimizable} = 0;

  return $self;
}

sub _new {
  my ( $class ) = @_;
  bless { data => [],
	  sum_x   => 0,
	  sum_x2  => 0,
	  sum_cnt => 0,
	  min => undef,
	  max => undef,
	  optimizable => 1 }, $class;
}

# =================================================================
# Accessors

sub count {
  my ( $self ) = @_;
  return $self->{sum_cnt};
}

sub range {
  my ( $self ) = @_;

  if( wantarray ) {
    return ( $self->{min}, $self->{max} );
  }

  return $self->{max};
}

sub extended_range {
  my ( $self ) = @_;

  my ( $min, $max ) = $self->range();
  my $w = $self->{extension}*$self->default_bandwidth();

  if( wantarray ) {
    return ( $min - $w, $max + $w );
  }

  return $max + $w;
}

sub default_bandwidth {
  my ( $self ) = @_;

  if( $self->{sum_cnt} == 0 ) { return undef; }

  my $x  = $self->{sum_x}/$self->{sum_cnt};
  my $x2 = $self->{sum_x2}/$self->{sum_cnt};
  my $sigma = sqrt( $x2 - $x**2 );

  # This is the optimal bandwidth if the point distribution is Gaussian.
  # (Applied Smoothing Techniques for Data Analysis
  # by Adrian W, Bowman & Adelchi Azzalini (1997)) */
  return $sigma * ( (3.0*$self->{sum_cnt}/4.0)**(-1.0/5.0) );
}

# =================================================================
# Adding Data

sub add_data {
  my ( $self, $x, $y, $w ) = @_;

  unless( _isNumber( $x ) ) { croak "Input ,$x, is not numeric."; }

  if( !defined( $y ) ) {
    $y = 1;
  } elsif( !_notNegative( $y ) ) {
    croak "Weight ,$y, must be non-negative number.";
  }

  if( defined( $w ) && !_isPositive( $w ) ) {
    croak "Bandwidth ,$w, must be strictly positive number in add_data.";
  }

  # If no bandwidth has been specified, $w will be undef!
  push @{ $self->{data} }, { pos => $x, cnt => $y, wid => $w };

  # Update summary statistics as we go along:
  $self->{sum_x}   += $y*$x;
  $self->{sum_x2}  += $y*$x*$x;
  $self->{sum_cnt} += $y;

  if( scalar @{ $self->{data} } == 1 ) {
    $self->{min} = $x;
    $self->{max} = $x;
  } else {
    $self->{min} = $x < $self->{min} ? $x : $self->{min};
    $self->{max} = $x > $self->{max} ? $x : $self->{max};
  }

  return;
}

# =================================================================
# Kernel Estimate

sub pdf {
  my $self = shift;
  return $self->_impl( 'pdf', 'default', @_ );
}

sub pdf_width_from_data {
  my $self = shift;
  return $self->_impl( 'pdf', 'fromdata', @_ );
}

# sub pdf_optimal {
#   my $self = shift;
#   return $self->_impl( 'pdf', 'optimal', @_ );
# }


sub cdf {
  my $self = shift;
  return $self->_impl( 'cdf', 'default', @_ );
}

sub cdf_width_from_data {
  my $self = shift;
  return $self->_impl( 'cdf', 'fromdata', @_ );
}

# sub cdf_optimal {
#   my $self = shift;
#   return $self->_impl( 'cdf', 'optimal', @_ );
# }

sub _curvature {
  my $self = shift;
  return $self->_impl( 'curvature', 'default', @_ );
}

sub _impl {
  my ( $self, $mode, $bandwidth_mode, $x, $w ) = @_;

  unless( $mode eq 'pdf' || $mode eq 'cdf'
	  || $mode eq 'curvature' ) { die "Illegal mode: ,$mode,"; }
  unless( _isNumber( $x ) ) { croak "Position ,$x, must be numeric."; }

  # If no data (or only data w/ weight zero), return immediately
  my $count = $self->count();
  if( $count == 0 ) { return 0; }

  # If bandwidth is from data, calculate result and return immediately
  if( $bandwidth_mode eq 'fromdata' ) {
    my $y = 0;
    for my $p ( @{ $self->{data} } ) {
      unless( defined $p->{wid} ) {
	carp "Undefined bandwidth in data at position " . $p->{pos}
	  . ". Using default bandwidth.";
	$w = $self->default_bandwidth();
      } else {
	$w = $p->{wid};
      }

      $y += $p->{cnt} * $self->{$mode}( $x, $p->{pos}, $w );
    }
    return $y/$count;
  }

  # ... otherwise, determine bandwidth
  if( $bandwidth_mode eq 'default' ) {
    if( !defined( $w ) ) {
      $w = $self->default_bandwidth();
    } elsif( !_notNegative( $w ) ) {
      croak "Bandwidth ,$w, must be strictly positive number.";
    }

# } elsif( $bandwidth_mode eq 'optimal' ) {
#   $w = $self->optimal_bandwidth();
  }

  # ... now use bandwidth from above to find result
  my $y = 0;
  for my $q ( @{ $self->{data} } ) {
    $y += $q->{cnt} * $self->{$mode}( $x, $q->{pos}, $w );
  }

  return $y/$count;
}

# =================================================================
# Classical Histograms

sub histogram {
  my ( $self, $bins ) = @_;

  unless( _isPositive( $bins ) && $bins==int($bins) ) {
    croak "Number of bins must be strictly positive integer.";
  }

  if( $self->count() == 0 ) { return []; }

  my ( $min, $max ) = $self->range();

  if( $bins == 1 ) {
    return [ { pos => 0.5*($max-$min), cnt => $self->count() } ];
  }

  my $w = ($max - $min)/($bins - 1);

  my @histo = ();
  for my $k ( 0..$bins-1 ) {
    push @histo, { pos => $min + $k*$w, cnt => 0 };
  }

  for my $p ( @{ $self->{data} } ) {
    my $i = int( ($p->{pos} - ( $min - 0.5*$w ) )/$w );
    my ( $lo, $hi ) = ( $min + ($i-0.5)*$w, $min + ($i+0.5)*$w );

    my $x = $p->{pos};

    if(     $x < $lo )             { $i -= 1; }
    elsif( $lo <= $x && $x < $hi ) { $i = $i; }
    elsif( $hi <= $x )             { $i += 1; }

    $histo[ $i ]->{cnt} += $p->{cnt};
  }

  return \@histo;
}

sub distribution_function {
  my ( $self ) = @_;

  my @sorted = sort { $a->{pos} <=> $b->{pos} } @{ $self->{data} };

  my @dist = ();
  my $cumul = 0;
  for my $p ( @sorted ) {
    $cumul += $p->{cnt};
    push @dist, { pos => $p->{pos}, cnt => $cumul };
  }

  return \@dist;
}

# =================================================================
# Input validation

# In general: undef evaluates to invalid input!

sub _isNumber {
  my ( $in ) = @_;
  if( defined( $in ) &&
      $in =~ /^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$/ ) { return 1; }
  return 0;
}

sub _isPositive {
  my ( $in ) = @_;
  if( _isNumber( $in ) && $in > 0 ) { return 1; }
  return 0;
}

sub _notNegative {
  my ( $in ) = @_;
  if( _isNumber( $in ) && $in >= 0 ) { return 1; }
  return 0;
}

# =================================================================
# Optimal Bandwidth

# Development history:
#
# Equation solver:
# 1) Straight iteration
# 2) Newton's method
# 3) Sekant method
# 4) Bisection
#
# Integration:
# 1) Numerical differentiation
# 2) Numerical differentiation, with equal step width as integration
# 3) Symbolic differentiation
# 4) Adaptive step size integration
# 5) Romberg integration (not implemented)

# This routine solves the equation encoded in _optimal_bandwidth_equation
# using the secant method.

sub optimal_bandwidth {
  my $self = shift;
  my $n    = @_ ? shift : 25;
  my $eps  = @_ ? shift : 1e-3;

  unless( $self->{optimizable} ) {
    croak "Bandwidth Optimization not available for this type of kernel.";
  }

  if( $self->{sum_cnt} == 0 ) { return undef; }

  my $x0 = $self->default_bandwidth();
  my $y0 = $self->_optimal_bandwidth_equation( $x0 );

  my $x = 0.8*$x0;
# my $x = $x0 * ( 1 - 1e-6 );
  my $y = $self->_optimal_bandwidth_equation( $x );

  my $dx = 0;

  my $i = 0;
  while( $i++ < $n ) {
    $x -= $y*($x0-$x)/($y0-$y);
    $y = $self->_optimal_bandwidth_equation( $x );

    if( abs($y) < $eps*$y0 ) { last }
  }

  if( wantarray ) { return ( $x, $i ); }
  return $x;
}

# This routine uses the secant method.

sub optimal_bandwidth_safe {
  my $self = shift;
  my $x0 = @_ ? shift : $self->default_bandwidth() / $self->count();
  my $x1 = @_ ? shift : 2*$self->default_bandwidth();
  my $eps  = @_ ? shift : 1e-3;

  unless( $self->{optimizable} ) {
    croak "Bandwidth Optimization not available for this type of kernel.";
  }

  if( $self->{sum_cnt} == 0 ) { return undef; }

  my $y0 = $self->_optimal_bandwidth_equation( $x0 );
  my $y1 = $self->_optimal_bandwidth_equation( $x1 );

  unless( $y0 * $y1 < 0 ) {
    croak "Interval [ f(x0=$x0)=$y0 : f(x1=$x1)=$y1 ] does not bracket root.";
  }

  my ( $x, $y, $i ) = ( 0, 0, 0 );
  while( abs( $x0 - $x1 ) > $eps*$x1 ) {
    $i += 1;

    $x = ( $x0 + $x1 )/2;
    $y = $self->_optimal_bandwidth_equation( $x );

    if( abs( $y ) < $eps*$y0 ) { last }

    if( $y * $y0 < 0 ) {
      ( $x1, $y1 ) = ( $x, $y );
    } else {
      ( $x0, $y0 ) = ( $x, $y );
    }
  }

  if( wantarray ) { return ( $x, $i ); }
  return $x;
}

# This routine encodes the self-consistent equation that is fulfilled
# by the optimal bandwidth. Notation according to Bowman & Azzalini.

sub _optimal_bandwidth_equation {
  my ( $self, $w ) = @_;

  my $alpha = 1.0/(2.0*sqrt( 3.14159265358979323846 ) );
  my $sigma = 1.0;
  my $n = $self->count();
  my $q = $self->_stiffness_integral( $w );

  return $w - ( ($n*$q*$sigma**4)/$alpha )**(-1.0/5.0);
}

# This routine calculates the integral over the square of the curvature
# (it: Int (f'')**2 ) using the trapezoidal rule. The routine decreases
# the step size by half until the relative error in the last step is less
# than epsilon.

sub _stiffness_integral {
  my ( $self, $w ) = @_;

  my $eps = 1e-4;

  my ( $mn, $mx ) = $self->extended_range();
  my $n = 1;
  my $dx = ($mx-$mn)/$n;

  my $yy = 0.5*($self->_curvature($mn,$w)**2+$self->_curvature($mx,$w)**2)*$dx;

  # The trapezoidal rule guarantees a relative error of better than eps
  # for some number of steps less than maxn.
  my $maxn = ($mx-$mn)/sqrt($eps);

  # This is not ideal, but I want to cap the total computation spent here:
  $maxn = ( $maxn > 2048 ? 2048 : $maxn );

  for( my $n=2; $n<=$maxn; $n*=2 ) {
    $dx /= 2.0;

    my $y = 0;
    for( my $i=1; $i<=$n-1; $i+=2 ) {
      $y += $self->_curvature( $mn + $i*$dx, $w )**2;
    }
    $yy = 0.5*$yy + $y*$dx;

    # Make at least 8 steps, then evaluate the relative change between steps
    if( $n > 8 && abs($y*$dx-0.5*$yy) < $eps*$yy ) { last }
  }

  return $yy;
}

# =================================================================
# Kernels

sub _gauss_pdf {
  my ( $x, $m, $s ) = @_;
  my $z = ($x - $m)/$s;
  return exp(-0.5*$z*$z)/( $s*sqrt( 2.0*3.14159265358979323846 ) );
}

# Abramowitz & Stegun, 26.2.17
sub _gauss_cdf {
  my ( $x, $m, $s ) = @_;

  my $z = abs( $x - $m)/$s;
  my $t = 1.0/(1.0 + 0.2316419*$z);
  my $y = $t*( 0.319381530
	       + $t*( -0.356563782
		    + $t*( 1.781477937
			  + $t*( -1.821255978 + $t*1.330274429 ) ) ) );

  if( $x >= $m ) {
    return 1.0 - _gauss_pdf( $x, $m, $s )*$y*$s;
  } else {
    return _gauss_pdf( $x, $m, $s )*$y*$s;
  }
}

sub _gauss_curvature {
  my ( $x, $m, $s ) = @_;
  my $z = ($x - $m)/$s;
  return ($z**2 - 1.0)*_gauss_pdf( $x, $m, $s )/$s**2;
}

sub _box_pdf {
  my ( $x, $m, $s ) = @_;
  if( $x < $m-0.5*$s || $x > $m+0.5*$s ) { return 0.0; }
  return 1.0/$s;
}

sub _box_cdf {
  my ( $x, $m, $s ) = @_;
  if( $x < $m-0.5*$s ) { return 0.0; }
  if( $x > $m+0.5*$s ) { return 1.0; }
  return ( $x-$m )/$s + 0.5;
}

sub _box_curvature {
  return 0;
}

sub _epanechnikov_pdf {
  my ( $x, $m, $s ) = @_;
  my $z = ($x-$m)/$s;
  if( abs($z) > 1 ) { return 0.0; }
  return 0.75*(1-$z**2)/$s;
}

sub _epanechnikov_cdf {
  my ( $x, $m, $s ) = @_;
  my $z = ($x-$m)/$s;
  if( $z < -1 ) { return 0.0; }
  if( $z >  1 ) { return 1.0; }
  return 0.25*(2.0 + 3.0*$z - $z**3 );
}

sub _epanechnikov_curvature {
  my ( $x, $m, $s ) = @_;
  my $z = ($x-$m)/$s;
  if( abs($z) > 1 ) { return 0; }
  return -1.5/$s**3;
}

1;

__END__