view release on metacpan or  search on metacpan

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

    $|++;
    my $self = shift;
    die "\nyou must first call read_file_for_bounding_box() before you can call\n" .   
        "metropolis_hastings()$!\n" unless $self->{_bounding_box};
    die "\nyou must first call normalize_input_histogram() before you can call\n" .
        "metropolis_hastings()$!\n" unless $self->{_normalized_input_hist};
    die "\nyou must first call set_sigmas_for_proposal_density() before you can call\n" .
        "metropolis_hastings()$!\n" unless $self->{_sigmax_for_proposal_density};
    my $box = $self->{_bounding_box};
    my $N_discard = $self->{_how_many_to_discard};
    my $N_iterations = $self->{_N} + $N_discard;
    my $sample = [$box->[0][0] + ($box->[0][1] - $box->[0][0]) / 2.0,
                  $box->[1][0] + ($box->[1][1] - $box->[1][0]) / 2.0];
    while ($self->desired_density($sample) == 0) {
        $sample = [$self->{_bounding_box}->[0][0] +
             random_uniform() * ($self->{_bounding_box}->[0][1] - $self->{_bounding_box}->[0][0]),
                   $self->{_bounding_box}->[1][0] +
             random_uniform() * ($self->{_bounding_box}->[1][1] - $self->{_bounding_box}->[1][0])]
    }
    print "\nstarting sample: @$sample\n" unless $self->{_command_line_mode};
    if ($self->{_command_line_mode}) {
        print "\nThe Metropolis-Hastings algorithm will be run over 2500 iterations.\n" .
              "Of the 2500 points generated, the first 500 will be discarded.\n" .
              "Each dot shown below stands for 50 iterations of the algorithm.\n\n";
    } else {
        print "\nThe Metropolis-Hastings algorithm will be run for $N_discard more iterations\n" .
              "than the number of points you requested.  The first $N_discard initial points thus\n" .
              "generated will be discarded from the final output.\n\n";
    }
    my @arr;
    foreach my $i (0..$N_iterations-1) {
        unless ($self->{_command_line_mode}) {
            print "\nIteration number: $i  (out of $N_iterations)\n" if $i % ($N_iterations / 10) == 0;
        }
        print ". " if $i % 50 == 0;
        # Get proposal probability q( $y | $x ).
        my ($newsample, $prob) = $self->get_sample_using_proposal( $sample ); 
        my $a1 = $self->desired_density( $newsample ) / $self->desired_density( $sample );
        my $a2 = $self->proposal_density( $sample, $newsample ) / $prob;
        my $a = $a1 * $a2;
        my $u = random_uniform();
        if ( $a >= 1 ) {
            $sample = $newsample;
            push @arr, $sample;
        } elsif ($u < $a) {
            $sample = $newsample;
            push @arr, $sample;
        } else {
            push @arr, $sample;
        }
    }
    print "\nTotal number of iterations run: $N_iterations\n\n" unless $self->{_command_line_mode};
    $self->{_generated_points} = \@arr;
}

# From the standpoint of matching the input histogram, the quality of the random
# points generated by the Metropolis-Hastings algorithm is affected by what sort of a
# probability distribution is used for the proposal density function.  This module
# uses a bivariate Gaussian density function for this purpose whose widths along x
# and y are controlled by constructor parameter "proposal_density_width" that
# defaults to 0.1.  The function here sets the values for the sigmax and sigmay
# parameters of this bivariate density.