#!/usr/local/bin/webwork-perl

=head1 Statistics Macros

=head3 Normal distribution

=pod 

	Usage: normal(a, b, mean=>0, deviation=>1); 

Computes the probability of x being in the interval (a,b) for normal distribution.
The first two arguments are required. Use '-infty' for negative infinity, and 'infty' or '+infty' for positive infinity. 
The mean and deviation are optional, and are 0 and 1 respectively by default.

=cut

sub normal { 
	my $a = shift;
	my $b = shift;
        my %options = @_;
        my $mean = $options{'mean'} if defined ($options{'mean'});
        $mean = 0 unless defined $mean;
        my $deviation = $options{'deviation'} if defined ($options{'deviation'});
	$deviation = 1 unless defined $deviation;
	if ($deviation <= 0) {
		warn 'Deviation must be a positive number.';
		return;
		}
	if ( $a eq '-infty' ) { $a = -6*$deviation + $mean; }
        if ( $b eq 'infty' ) { $b = 6*$deviation + $mean; }
	if ( $b eq '+infty' ) { $b = 6*$deviation + $mean; }
        my $z_score_of_a = ($a - $mean)/$deviation;
        my $z_score_of_b = ($b - $mean)/$deviation;
        my $function = sub { my $x=shift;
                             $E**(-$x**2/2)/sqrt(2*$PI);
                             }; 
        my $prob = romberg($function, $z_score_of_a, $z_score_of_b);
        $prob;
}

=head3 "Inverse" of normal distribution

=pod 

	Usage: inv_normal(prob, mean=>0, deviation=>1);

Computes the positive number b such that the probability of x being in the interval (0,b)
is equal to the given probability (first argument). The mean and deviation are
optional, and are 0 and 1 respectively by default. 
Caution: since students may use tables, they may only be able to provide the answer correct to 2 or 3
decimal places. Use tolerance when evaluating answers.

=cut

sub inv_normal { 
	my $prob = shift;
	my %options = @_;
        my $mean = $options{'mean'} if defined ($options{'mean'});
        $mean = 0 unless defined $mean;
        my $deviation = $options{'deviation'} if defined ($options{'deviation'});
        $deviation = 1 unless defined $deviation;
        if ($deviation <= 0) {
                warn 'Deviation must be a positive number.';
                return;
                }
        my $function = sub { my $x=shift;
                             $E**(-$x**2/2)/sqrt(2*$PI);
                             };
        my $z_score_of_b = inv_romberg($function, 0, $prob);
        my $b = $z_score_of_b * $deviation + $mean;
        $b;
}    

##########################################

1;