view release on metacpan or  search on metacpan

lib/Algorithm/CurveFit.pm  view on Meta::CPAN

package Algorithm::CurveFit;

use 5.006;
use strict;
use warnings;

our $VERSION = '1.06';

require Exporter;

our @ISA = qw(Exporter);

our %EXPORT_TAGS = (
    'all' => [ qw( curve_fit ) ]
);

our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );

our @EXPORT = qw();

use Carp qw/confess/;
use Math::Symbolic qw/parse_from_string/;
use Math::MatrixReal;
use Data::Dumper;

# machine epsilon
use constant EPS => 2.2e-16;
use constant SQRT_EPS => sqrt(EPS);

sub curve_fit {
    shift @_ if not ref $_[0] and defined $_[0] and $_[0] eq 'Algorithm::CurveFit';

    confess('Uneven number of arguments to Algorithm::CurveFit::curve_fit.')
      if @_ % 2;

    my %args = @_;

    # Formula
    confess("Missing 'formula' parameter.") if not defined $args{formula};
    my $formula;
    if (ref($args{formula}) =~ /^Math::Symbolic/) {
        $formula = $args{formula};
    }
    else {
        eval { $formula = parse_from_string( $args{formula} ); };
        confess( "Cannot parse formula '" . $args{formula} . "'. ($@)" )
          if not defined $formula or $@;
    }

    # Variable (optional)
    my $variable = $args{variable};
    $variable = 'x' if not defined $variable;
    confess("Formula '"
          . $args{formula}
          . "' not explicitly dependent on "
          . "variable '$variable'." )
      if not grep { $_ eq $variable } $formula->explicit_signature();

    # Parameters
    my $params = $args{params};

    confess("Parameter 'params' has to be an array reference.")
      if not defined $params
      or not ref($params) eq 'ARRAY';
    my @parameters = @$params;
    confess('No parameters specified.') if not @parameters;
    confess('Individual parameters need to be array references.')
      if grep { not defined $_ or not ref($_) eq 'ARRAY' } @parameters;
    foreach my $p (@parameters) {
        confess("Weird parameter\n'"
              . Dumper($p)
              . "' Should have the format\n"
              . "[ NAME_STRING, GUESSED_VALUE, ACCURACY ]\n"
              . "With the accuracy being optional. See docs." )
          if @$p > 3
          or @$p < 2
          or grep { not defined $_ } @$p;

        confess("Formula '"
              . $args{formula}
              . "' not explicitly dependent on "
              . "parameter '"
              . $p->[0]
              . "'." )

lib/Algorithm/CurveFit.pm  view on Meta::CPAN

        # solve
        my $LR = $matrix->decompose_LR();
        my ( $dim, $x, $B ) = $LR->solve_LR($vector);

        # extract parameter modifications and test for convergence
        my $last = 1;
        foreach my $pno ( 1 .. @parameters ) {
            my $dlambda = $x->element( $pno, 1 );
            $last = 0 if abs($dlambda) > $parameters[ $pno - 1 ][2];
            $parameters[ $pno - 1 ][1] += $dlambda;
        }
        last if $last;
    }

    # Recalculate dbeta for the squared residuals.
    my @beta =
      map {
        $ydata[$_] - $formula_sub->(
            $xdata[$_],
            map { $_->[1] } @parameters
          )
      } 0 .. $#xdata;
    $dbeta = Math::MatrixReal->new_from_cols( [ \@beta ] );

    my $square_residual = $dbeta->scalar_product($dbeta);
    return $square_residual;
}

1;
__END__

=head1 NAME

Algorithm::CurveFit - Nonlinear Least Squares Fitting

=head1 SYNOPSIS

  use Algorithm::CurveFit;

  # Known form of the formula
  my $formula = 'c + a * x^2';
  my $variable = 'x';
  my @xdata = read_file('xdata'); # The data corresponsing to $variable
  my @ydata = read_file('ydata'); # The data on the other axis
  my @parameters = (
      # Name    Guess   Accuracy
      ['a',     0.9,    0.00001],  # If an iteration introduces smaller
      ['c',     20,     0.00005],  # changes that the accuracy, end.
  );
  my $max_iter = 100; # maximum iterations

  my $square_residual = Algorithm::CurveFit->curve_fit(
      formula            => $formula, # may be a Math::Symbolic tree instead
      params             => \@parameters,
      variable           => $variable,
      xdata              => \@xdata,
      ydata              => \@ydata,
      maximum_iterations => $max_iter,
  );

  use Data::Dumper;
  print Dumper \@parameters;
  # Prints
  # $VAR1 = [
  #          [
  #            'a',
  #            '0.201366784209602',
  #            '1e-05'
  #          ],
  #          [
  #            'c',
  #            '1.94690440147554',
  #            '5e-05'
  #          ]
  #        ];
  #
  # Real values of the parameters (as demonstrated by noisy input data):
  # a = 0.2
  # c = 2

=head1 DESCRIPTION

C<Algorithm::CurveFit> implements a nonlinear least squares curve fitting
algorithm. That means, it fits a curve of known form (sine-like, exponential,
polynomial of degree n, etc.) to a given set of data points.

For details about the algorithm and its capabilities and flaws, you're
encouraged to read the MathWorld page referenced below. Note, however, that it
is an iterative algorithm that improves the fit with each iteration until it
converges. The following rule of thumb usually holds true:

=over 2

=item

A good guess improves the probability of convergence and the quality
of the fit.

=item

Increasing the number of free parameters decreases the quality and
convergence speed.

=item

Make sure that there are no correlated parameters such as in 'a + b * e^(c+x)'.
(The example can be rewritten as 'a + b * e^c * e^x' in which 'c' and 'b' are
basically equivalent parameters.

=back

The curve fitting algorithm is accessed via the 'curve_fit' subroutine.
It requires the following parameters as 'key => value' pairs:

=over 2

=item formula

The formula should be a string that can be parsed by Math::Symbolic.
Alternatively, it can be an existing Math::Symbolic tree.
Please refer to the documentation of that module for the syntax.