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
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
|
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016 The GSL Team.
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with the
Invariant Sections being "GNU General Public License" and "Free Software
Needs Free Documentation", the Front-Cover text being "A GNU Manual",
and with the Back-Cover Text being (a) (see below). A copy of the
license is included in the section entitled "GNU Free Documentation
License".
(a) The Back-Cover Text is: "You have the freedom to copy and modify this
GNU Manual." -->
<!-- Created by GNU Texinfo 5.1, http://www.gnu.org/software/texinfo/ -->
<head>
<title>GNU Scientific Library – Reference Manual: VEGAS</title>
<meta name="description" content="GNU Scientific Library – Reference Manual: VEGAS">
<meta name="keywords" content="GNU Scientific Library – Reference Manual: VEGAS">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<link href="index.html#Top" rel="start" title="Top">
<link href="Function-Index.html#Function-Index" rel="index" title="Function Index">
<link href="Monte-Carlo-Integration.html#Monte-Carlo-Integration" rel="up" title="Monte Carlo Integration">
<link href="Monte-Carlo-Examples.html#Monte-Carlo-Examples" rel="next" title="Monte Carlo Examples">
<link href="MISER.html#MISER" rel="previous" title="MISER">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.indentedblock {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smallindentedblock {margin-left: 3.2em; font-size: smaller}
div.smalllisp {margin-left: 3.2em}
kbd {font-style:oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nocodebreak {white-space:nowrap}
span.nolinebreak {white-space:nowrap}
span.roman {font-family:serif; font-weight:normal}
span.sansserif {font-family:sans-serif; font-weight:normal}
ul.no-bullet {list-style: none}
-->
</style>
</head>
<body lang="en" bgcolor="#FFFFFF" text="#000000" link="#0000FF" vlink="#800080" alink="#FF0000">
<a name="VEGAS"></a>
<div class="header">
<p>
Next: <a href="Monte-Carlo-Examples.html#Monte-Carlo-Examples" accesskey="n" rel="next">Monte Carlo Examples</a>, Previous: <a href="MISER.html#MISER" accesskey="p" rel="previous">MISER</a>, Up: <a href="Monte-Carlo-Integration.html#Monte-Carlo-Integration" accesskey="u" rel="up">Monte Carlo Integration</a> [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="VEGAS-1"></a>
<h3 class="section">25.4 VEGAS</h3>
<a name="index-VEGAS-Monte-Carlo-integration"></a>
<a name="index-importance-sampling_002c-VEGAS"></a>
<p>The <small>VEGAS</small> algorithm of Lepage is based on importance sampling. It
samples points from the probability distribution described by the
function <em>|f|</em>, so that the points are concentrated in the regions
that make the largest contribution to the integral.
</p>
<p>In general, if the Monte Carlo integral of <em>f</em> is sampled with
points distributed according to a probability distribution described by
the function <em>g</em>, we obtain an estimate <em>E_g(f; N)</em>,
</p>
<div class="example">
<pre class="example">E_g(f; N) = E(f/g; N)
</pre></div>
<p>with a corresponding variance,
</p>
<div class="example">
<pre class="example">\Var_g(f; N) = \Var(f/g; N).
</pre></div>
<p>If the probability distribution is chosen as <em>g = |f|/I(|f|)</em> then
it can be shown that the variance <em>V_g(f; N)</em> vanishes, and the
error in the estimate will be zero. In practice it is not possible to
sample from the exact distribution <em>g</em> for an arbitrary function, so
importance sampling algorithms aim to produce efficient approximations
to the desired distribution.
</p>
<p>The <small>VEGAS</small> algorithm approximates the exact distribution by making a
number of passes over the integration region while histogramming the
function <em>f</em>. Each histogram is used to define a sampling
distribution for the next pass. Asymptotically this procedure converges
to the desired distribution. In order
to avoid the number of histogram bins growing like <em>K^d</em> the
probability distribution is approximated by a separable function:
<em>g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ...</em>
so that the number of bins required is only <em>Kd</em>.
This is equivalent to locating the peaks of the function from the
projections of the integrand onto the coordinate axes. The efficiency
of <small>VEGAS</small> depends on the validity of this assumption. It is most
efficient when the peaks of the integrand are well-localized. If an
integrand can be rewritten in a form which is approximately separable
this will increase the efficiency of integration with <small>VEGAS</small>.
</p>
<p><small>VEGAS</small> incorporates a number of additional features, and combines both
stratified sampling and importance sampling. The integration region is
divided into a number of “boxes”, with each box getting a fixed
number of points (the goal is 2). Each box can then have a fractional
number of bins, but if the ratio of bins-per-box is less than two, Vegas switches to a
kind variance reduction (rather than importance sampling).
</p>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005falloc"></a>Function: <em>gsl_monte_vegas_state *</em> <strong>gsl_monte_vegas_alloc</strong> <em>(size_t <var>dim</var>)</em></dt>
<dd><a name="index-gsl_005fmonte_005fvegas_005fstate"></a>
<p>This function allocates and initializes a workspace for Monte Carlo
integration in <var>dim</var> dimensions. The workspace is used to maintain
the state of the integration.
</p></dd></dl>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005finit"></a>Function: <em>int</em> <strong>gsl_monte_vegas_init</strong> <em>(gsl_monte_vegas_state* <var>s</var>)</em></dt>
<dd><p>This function initializes a previously allocated integration state.
This allows an existing workspace to be reused for different
integrations.
</p></dd></dl>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005fintegrate"></a>Function: <em>int</em> <strong>gsl_monte_vegas_integrate</strong> <em>(gsl_monte_function * <var>f</var>, double <var>xl</var>[], double <var>xu</var>[], size_t <var>dim</var>, size_t <var>calls</var>, gsl_rng * <var>r</var>, gsl_monte_vegas_state * <var>s</var>, double * <var>result</var>, double * <var>abserr</var>)</em></dt>
<dd><p>This routines uses the <small>VEGAS</small> Monte Carlo algorithm to integrate the
function <var>f</var> over the <var>dim</var>-dimensional hypercubic region
defined by the lower and upper limits in the arrays <var>xl</var> and
<var>xu</var>, each of size <var>dim</var>. The integration uses a fixed number
of function calls <var>calls</var>, and obtains random sampling points using
the random number generator <var>r</var>. A previously allocated workspace
<var>s</var> must be supplied. The result of the integration is returned in
<var>result</var>, with an estimated absolute error <var>abserr</var>. The result
and its error estimate are based on a weighted average of independent
samples. The chi-squared per degree of freedom for the weighted average
is returned via the state struct component, <var>s->chisq</var>, and must be
consistent with 1 for the weighted average to be reliable.
</p></dd></dl>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005ffree"></a>Function: <em>void</em> <strong>gsl_monte_vegas_free</strong> <em>(gsl_monte_vegas_state * <var>s</var>)</em></dt>
<dd><p>This function frees the memory associated with the integrator state
<var>s</var>.
</p></dd></dl>
<p>The <small>VEGAS</small> algorithm computes a number of independent estimates of the
integral internally, according to the <code>iterations</code> parameter
described below, and returns their weighted average. Random sampling of
the integrand can occasionally produce an estimate where the error is
zero, particularly if the function is constant in some regions. An
estimate with zero error causes the weighted average to break down and
must be handled separately. In the original Fortran implementations of
<small>VEGAS</small> the error estimate is made non-zero by substituting a small
value (typically <code>1e-30</code>). The implementation in GSL differs from
this and avoids the use of an arbitrary constant—it either assigns
the value a weight which is the average weight of the preceding
estimates or discards it according to the following procedure,
</p>
<dl compact="compact">
<dt>current estimate has zero error, weighted average has finite error</dt>
<dd>
<p>The current estimate is assigned a weight which is the average weight of
the preceding estimates.
</p>
</dd>
<dt>current estimate has finite error, previous estimates had zero error</dt>
<dd>
<p>The previous estimates are discarded and the weighted averaging
procedure begins with the current estimate.
</p>
</dd>
<dt>current estimate has zero error, previous estimates had zero error</dt>
<dd>
<p>The estimates are averaged using the arithmetic mean, but no error is computed.
</p></dd>
</dl>
<p>The convergence of the algorithm can be tested using the overall
chi-squared value of the results, which is available from the
following function:
</p>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005fchisq"></a>Function: <em>double</em> <strong>gsl_monte_vegas_chisq</strong> <em>(const gsl_monte_vegas_state * <var>s</var>)</em></dt>
<dd><p>This function returns the chi-squared per degree of freedom for the
weighted estimate of the integral. The returned value should be close
to 1. A value which differs significantly from 1 indicates that the
values from different iterations are inconsistent. In this case the
weighted error will be under-estimated, and further iterations of the
algorithm are needed to obtain reliable results.
</p></dd></dl>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005frunval"></a>Function: <em>void</em> <strong>gsl_monte_vegas_runval</strong> <em>(const gsl_monte_vegas_state * <var>s</var>, double * <var>result</var>, double * <var>sigma</var>)</em></dt>
<dd><p>This function returns the raw (unaveraged) values of the integral
<var>result</var> and its error <var>sigma</var> from the most recent iteration
of the algorithm.
</p></dd></dl>
<p>The <small>VEGAS</small> algorithm is highly configurable. Several parameters
can be changed using the following two functions.
</p>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005fparams_005fget"></a>Function: <em>void</em> <strong>gsl_monte_vegas_params_get</strong> <em>(const gsl_monte_vegas_state * <var>s</var>, gsl_monte_vegas_params * <var>params</var>) </em></dt>
<dd><p>This function copies the parameters of the integrator state into the
user-supplied <var>params</var> structure.
</p></dd></dl>
<dl>
<dt><a name="index-gsl_005fmonte_005fvegas_005fparams_005fset"></a>Function: <em>void</em> <strong>gsl_monte_vegas_params_set</strong> <em>(gsl_monte_vegas_state * <var>s</var>, const gsl_monte_vegas_params * <var>params</var>) </em></dt>
<dd><p>This function sets the integrator parameters based on values provided
in the <var>params</var> structure.
</p></dd></dl>
<p>Typically the values of the parameters are first read using
<code>gsl_monte_vegas_params_get</code>, the necessary changes are made to
the fields of the <var>params</var> structure, and the values are copied
back into the integrator state using
<code>gsl_monte_vegas_params_set</code>. The functions use the
<code>gsl_monte_vegas_params</code> structure which contains the following
fields:
</p>
<dl>
<dt><a name="index-alpha-1"></a>Variable: <em>double</em> <strong>alpha</strong></dt>
<dd><p>The parameter <code>alpha</code> controls the stiffness of the rebinning
algorithm. It is typically set between one and two. A value of zero
prevents rebinning of the grid. The default value is 1.5.
</p></dd></dl>
<dl>
<dt><a name="index-iterations"></a>Variable: <em>size_t</em> <strong>iterations</strong></dt>
<dd><p>The number of iterations to perform for each call to the routine. The
default value is 5 iterations.
</p></dd></dl>
<dl>
<dt><a name="index-stage"></a>Variable: <em>int</em> <strong>stage</strong></dt>
<dd><p>Setting this determines the <em>stage</em> of the calculation. Normally,
<code>stage = 0</code> which begins with a new uniform grid and empty weighted
average. Calling <small>VEGAS</small> with <code>stage = 1</code> retains the grid from the
previous run but discards the weighted average, so that one can “tune”
the grid using a relatively small number of points and then do a large
run with <code>stage = 1</code> on the optimized grid. Setting <code>stage =
2</code> keeps the grid and the weighted average from the previous run, but
may increase (or decrease) the number of histogram bins in the grid
depending on the number of calls available. Choosing <code>stage = 3</code>
enters at the main loop, so that nothing is changed, and is equivalent
to performing additional iterations in a previous call.
</p></dd></dl>
<dl>
<dt><a name="index-mode"></a>Variable: <em>int</em> <strong>mode</strong></dt>
<dd><p>The possible choices are <code>GSL_VEGAS_MODE_IMPORTANCE</code>,
<code>GSL_VEGAS_MODE_STRATIFIED</code>, <code>GSL_VEGAS_MODE_IMPORTANCE_ONLY</code>.
This determines whether <small>VEGAS</small> will use importance sampling or
stratified sampling, or whether it can pick on its own. In low
dimensions <small>VEGAS</small> uses strict stratified sampling (more precisely,
stratified sampling is chosen if there are fewer than 2 bins per box).
</p></dd></dl>
<dl>
<dt><a name="index-verbose"></a>Variable: <em>int</em> <strong>verbose</strong></dt>
<dt><a name="index-ostream"></a>Variable: <em>FILE *</em> <strong>ostream</strong></dt>
<dd><p>These parameters set the level of information printed by <small>VEGAS</small>. All
information is written to the stream <var>ostream</var>. The default setting
of <var>verbose</var> is <code>-1</code>, which turns off all output. A
<var>verbose</var> value of <code>0</code> prints summary information about the
weighted average and final result, while a value of <code>1</code> also
displays the grid coordinates. A value of <code>2</code> prints information
from the rebinning procedure for each iteration.
</p></dd></dl>
<p>The above fields and the <var>chisq</var> value can also be accessed
directly in the <code>gsl_monte_vegas_state</code> but such use is
deprecated.
</p>
<hr>
<div class="header">
<p>
Next: <a href="Monte-Carlo-Examples.html#Monte-Carlo-Examples" accesskey="n" rel="next">Monte Carlo Examples</a>, Previous: <a href="MISER.html#MISER" accesskey="p" rel="previous">MISER</a>, Up: <a href="Monte-Carlo-Integration.html#Monte-Carlo-Integration" accesskey="u" rel="up">Monte Carlo Integration</a> [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
</body>
</html>
|
