1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
use strict;
use warnings;
our $VERSION = '2.096';
pp_addpm({At=>'Top'},<<'EOD');
use strict;
use warnings;
=head1 NAME
PDL::GSL::CDF - PDL interface to GSL Cumulative Distribution Functions
=head1 DESCRIPTION
This is an interface to the Cumulative Distribution Function package present in the GNU Scientific Library.
Let us have a continuous random number distributions are defined by a probability density function C<p(x)>.
The cumulative distribution function for the lower tail C<P(x)> is defined by the integral of C<p(x)>, and
gives the probability of a variate taking a value less than C<x>. These functions are named B<cdf_NNNNNNN_P()>.
The cumulative distribution function for the upper tail C<Q(x)> is defined by the integral of C<p(x)>, and
gives the probability of a variate taking a value greater than C<x>. These functions are named B<cdf_NNNNNNN_Q()>.
The upper and lower cumulative distribution functions are related by C<P(x) + Q(x) = 1> and
satisfy C<0 E<lt>= P(x) E<lt>= 1> and C<0 E<lt>= Q(x) E<lt>= 1>.
The inverse cumulative distributions, C<x = Pinv(P)> and C<x = Qinv(Q)> give the values of C<x> which correspond
to a specific value of C<P> or C<Q>. They can be used to find confidence limits from probability values.
These functions are named B<cdf_NNNNNNN_Pinv()> and B<cdf_NNNNNNN_Qinv()>.
For discrete distributions the probability of sampling the integer value C<k> is given by C<p(k)>, where
C<sum_k p(k) = 1>. The cumulative distribution for the lower tail C<P(k)> of a discrete distribution is
defined as, where the sum is over the allowed range of the distribution less than or equal to C<k>.
The cumulative distribution for the upper tail of a discrete distribution C<Q(k)> is defined as giving the sum
of probabilities for all values greater than C<k>. These two definitions satisfy the identity C<P(k) + Q(k) = 1>.
If the range of the distribution is C<1> to C<n> inclusive then C<P(n) = 1>, C<Q(n) = 0>
while C<P(1) = p(1)>, C<Q(1) = 1 - p(1)>.
=head1 SYNOPSIS
use PDL;
use PDL::GSL::CDF;
my $p = gsl_cdf_tdist_P( $t, $df );
my $t = gsl_cdf_tdist_Pinv( $p, $df );
=cut
EOD
pp_addhdr('
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
static char *funcname;
static void cdf_error_handler(const char *reason, const char *file, int line, int status) {
char buf[200];
sprintf(buf,"Error in %s: %s", funcname, gsl_strerror(status));
barf(buf);
}
');
use Text::ParseWords qw(shellwords);
chomp(my $header = `gsl-config --cflags`);
$header =~ s#\\#/#g; # win32
($header) = map {s/^-I//;$_} grep /^-I/, shellwords $header;
$header .= '/gsl/gsl_cdf.h';
my $h = do { open my $fh, "<", $header or die "$header: $!"; local $/ = undef; <$fh> };
my @functions = $h =~ m/(double\s+gsl_cdf_.+?\)\s*;\s*)\n/xmsg;
my @func_defs;
my %p_type = (
'double' => 'double',
'unsigned int' => 'ulonglong',
);
my %func_desc = (
gsl_cdf_gaussian => ['The Gaussian Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Gaussian distribution with standard deviation I<sigma>.'],
gsl_cdf_ugaussian => ['The Unit Gaussian Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the unit Gaussian distribution.'],
gsl_cdf_exponential => ['The Exponential Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the exponential distribution with mean I<mu>.'],
gsl_cdf_laplace => ['The Laplace Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Laplace distribution with width I<a>.'],
gsl_cdf_exppow => ['The Exponential Power Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) for the exponential power distribution with parameters I<a> and I<b>.'],
gsl_cdf_cauchy => ['The Cauchy Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Cauchy distribution with scale parameter I<a>.'],
gsl_cdf_rayleigh => ['The Rayleigh Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Rayleigh distribution with scale parameter I<sigma>.'],
gsl_cdf_gamma => ['The Gamma Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the gamma distribution with parameters I<a> and I<b>.'],
gsl_cdf_flat => ['The Flat (Uniform) Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for a uniform distribution from I<a> to I<b>.'],
gsl_cdf_lognormal => ['The Lognormal Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the lognormal distribution with parameters I<zeta> and I<sigma>.'],
gsl_cdf_chisq => ['The Chi-squared Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the chi-squared distribution with I<nu> degrees of freedom.'],
gsl_cdf_fdist => ['The F-distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the F-distribution with I<nu1> and I<nu2> degrees of freedom.'],
gsl_cdf_tdist => ['The t-distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the t-distribution with I<nu> degrees of freedom.'],
gsl_cdf_beta => ['The Beta Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the beta distribution with parameters I<a> and I<b>.'],
gsl_cdf_logistic => ['The Logistic Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the logistic distribution with scale parameter I<a>.'],
gsl_cdf_pareto => ['The Pareto Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Pareto distribution with exponent I<a> and scale I<b>.'],
gsl_cdf_weibull => ['The Weibull Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Weibull distribution with scale I<a> and exponent I<b>.'],
gsl_cdf_gumbel1 => ['The Type-1 Gumbel Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-1 Gumbel distribution with parameters I<a> and I<b>.'],
gsl_cdf_gumbel2 => ['The Type-2 Gumbel Distribution', 'These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-2 Gumbel distribution with parameters I<a> and I<b>.'],
gsl_cdf_poisson => ['The Poisson Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the Poisson distribution with parameter I<mu>.'],
gsl_cdf_binomial => ['The Binomial Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the binomial distribution with parameters I<p> and I<n>.'],
gsl_cdf_negative_binomial => ['The Negative Binomial Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the negative binomial distribution with parameters I<p> and I<n>.'],
gsl_cdf_pascal => ['The Pascal Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the Pascal distribution with parameters I<p> and I<n>.'],
gsl_cdf_geometric => ['The Geometric Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the geometric distribution with parameter I<p>.'],
gsl_cdf_hypergeometric => ['The Hypergeometric Distribution', 'These functions compute the cumulative distribution functions P(k), Q(k) for the hypergeometric distribution with parameters I<n1>, I<n2> and I<t>.'],
);
for (@functions) {
s/\n\s+/ /xmsg;
s/const //g;
if (m/^(\w+)\s+(\w+)\s*\(\s*(.+)\s*\)\s*\;$/s) {
my ($out_type, $function, $pars) = ($1, $2, $3);
my @pars = split /,/, $pars;
for (@pars) {
if (m/^(.+)( \w+)$/) {
my ($type, $par) = ($1, $2);
s/^ | $//g for ($type, $par);
$par = lc $par;
$_ = [$p_type{$type}, $par];
}
}
push @func_defs, [ $out_type, $function, \@pars ];
}
}
@func_defs = sort { $a->[1] cmp $b->[1] } @func_defs; # sort by function name
for my $f (@func_defs) {
my ($out_type, $function, $pars) = @$f;
my $func_short = join '_', (split '_', $function)[0..2];
my ($p, $code) = print_ppdef($out_type, $function, @$pars);
my $desc = delete $func_desc{$func_short};
pp_addpm(qq{\n=head2 $desc->[0] (${func_short}_*)\n\n$desc->[1]\n\n=cut\n\n}) if $desc;
pp_def($function,
HandleBad => 1,
GenericTypes => ['D'],
Pars => $p,
Code => $code,
Doc => '',
);
}
sub print_ppdef {
my ($out_type, $function, @pars) = @_;
@pars = map $_->[0] eq 'double' ? [$_->[1]] : $_, @pars;
my $pars = join ' ', map "@$_();", @pars, ['[o]out'];
my $code = pp_line_numbers(__LINE__, <<EOF);
gsl_error_handler_t *current_handler;
funcname = "$function";
current_handler = gsl_set_error_handler(cdf_error_handler);
\$out() = $function(@{[join ', ', map "\$$_->[-1]()", @pars]});
gsl_set_error_handler(current_handler);
EOF
$code =
"PDL_IF_BAD(if ( !(" . join(' && ', map { "\$ISGOOD($_->[-1]())" } @pars ) . ") ) {" .
"\$SETBAD(out()); }" .
'else,) {' . $code . "}\n";
return ($pars, $code);
}
pp_addpm({At=>'Bot'},<<'EOD');
=head1 AUTHOR
Copyright (C) 2009 Maggie J. Xiong <maggiexyz users.sourceforge.net>
The GSL CDF module was written by J. Stover.
All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDL distribution.
=cut
EOD
pp_done();
|
