vstatmp_engl

other polygon in the circle and considering only the area between the polygons.

b) Importance Sampling

If there exists a majorant m(x) for the function y(x) to be integrated,

with the property that the indefinite integral M(x)

Z x

M(x) = m(x′)dx′

xa

can be inverted, we generate N0 x-values according to the distribution m(x). For each xi a further random number yi in the interval 0 < y < m(xi) is generated. Again, as for the simulation of distributions, points lying above y(xi) are rejected. The number N of the remaining events provides the integral

5.3.3 Weighting Method

a) Simple Weighting

We generate N random numbers xi in the interval xa < x < xb and average over the

function values:

XN

y = y(xi)/N .

i=1

An estimate for the integral (5.4) is given by

This method corresponds to the usual numerical integration, with the peculiarity that the supporting points on the abscissa are not chosen regularly but are distributed at random. This alone cannot be an advantage, and indeed the Monte Carlo integration in one and two dimensions for a given number of supporting points is less e cient than conventional methods. It is, however, superior to other methods for multi-dimensional integrations. Already in three dimensions it competes favorably in many cases.

To estimate the accuracy, we apply the usual statistical error estimation. We consider the numbers yi = y(xi) as N stochastic measurements of y. The expected mean squared error of y is then given by (4.3):

128 5 Monte Carlo Simulation

Fig. 5.13. Monte Carlo integration of the di erence between the function to be integrated and an integrable function.

The numerator is an estimate of the variance of the y distribution. The accuracy is the better, the smaller the fluctuations of the function around its mean value are.

b) Subtraction method

The accuracy can be improved through a reduction of the fluctuations of the integrand.

If we find a function ye(x) which is integrable analytically and does not di er too much from the original integrand y(x) we cast the integral into the form

We now have to evaluate by Monte Carlo only the second term with relatively small fluctuations (Fig. 5.13).

5.3.4 Reduction to Expected Values

In many cases it makes sense to factorize the integrand y(x) = f(x)y1(x) into a factor f(x) corresponding to a p.d.f. normalized inside the integration interval which

is easy to generate, and a second factor y1(x). To be e ective, the method requires that f is close to y. Our integral has now the form of an expected value:

= hy1i .

We generate values xi distributed according to f(x) and obtain from these an estimate for the integral I:

The estimate is again the better, the less the y1-values are fluctuating, i.e. the more similar the functions y and f are. The error estimate is analogous to (5.5).

5.3.5 Stratified Sampling

In stratified sampling the domain of integration is partitioned into sub-domains. Over each of these we integrate separately. The advantage is that the distribution in each sub-domain is more uniform and thus the fluctuations of the random variables are smaller and the statistical error is reduced. This method is somewhat antithetical to the basic idea of the simple Monte Carlo method, since it produces a more uniform (equidistant) distribution of the supporting points and requires some e ort to combine the errors from the di erent contributions. Thus we recommend it only if the integrand shows very strong variations.

5.4 General Remarks

Often we need to solve integrals over di erent domains but always with the same integrand. In these cases the Monte Carlo approach is particularly advantageous. We store all single simulated values (usually called “events”) and are able to select events afterwards according to the chosen domain, and obtain the integral with relatively small computing expense by summation. Similarly a change of event weights is possible without repeating the generation of the events.

Let us illustrate this feature with a mechanical example: If, for instance, we want to obtain the tensor of inertia for a complex mass distribution like a car, we distribute points stochastically within the body and store their coordinates together with the respective mass densities. With these data it is easy to calculate by summations the mass, the center of mass and the moments of inertia with respect to arbitrary axes. If desired, parts of the body can be eliminated simply by rejecting the corresponding points in the sums and di erent materials can be considered by changing the density.

In thermodynamic systems we are often interested in several mean values, like the mean free path length, mean kinetic or potential energy, velocities etc.. Once a

132 6 Parameter Inference I

but the likelihood alone is not su cient for a decision in favor of one or the other hypothesis. Additional considerations may play an important role. When the cards are located in a Swiss youth hostel we would consider the hypothesis Skat more sceptically than when the cards are found in a pub at Hamburg. We therefore would weight our hypotheses with prior probabilities (short: priors) which quantify this additional piece of information. Prior probabilities are often hard to estimate, often they are completely unknown. As a consequence, results depending on priors are model dependent.

We usually will avoid to introduce prior probabilities and stay with likelihoods but sometimes this is not possible. Then the results have to be interpreted conditional to the validity of the applied prior probabilities.

Similar situations as described in our trivial example also occur in empirical sciences. Whenever an observation is more or less compatible with two alternative theories, we cannot simply derive probabilities for the validities of theories based solely on the experimental data. Other criteria like the attractiveness of the theory, the compatibility with previous measurements will enter into our judgement. These additional attributes again can be taken into account by introducing prior probabilities which of course will depend to some extend on subjective prejudices. The cases where well founded quantitative priors are available are rare.

Some years ago, in an experiment a significant deviation from the gravitational law had been observed. From the two alternatives: H1: Newton’s law is correct and

H2: A 1/r2 term has to be added to the 1/r potential, the latter was much more compatible with the experimental data. In spite of this experimental evidence, hardly any serious physicist doubted the validity of the classical gravitational law. The reason is that our experience shows that the laws of nature are basically simple. In the mean time, as a consequence of the availability of more precise data, the 1/r2 hypothesis has been rejected.

Nevertheless it is correct and necessary to publish the observed data without weighting them with a prior, i.e. to restrict the presentation to the purely statistical result and to leave the subjective part to the reader of the publication.

The scientific method requires the possibility to compare and combine results from independent measurements. This is impossible when di erent authors apply varying priors. We will see that it is almost always possible to avoid the introduction of this kind of subjective assumptions into the statistical analysis.

In our second example the situation is even more complex because we are confronted with only one hypothesis and no well specified alternative. The validity of the alternative, e.g. a deviation from the equality of the distribution of the sexes is hardly measurable since an arbitrarily small deviation from the equality is present in any case. There is no other possibility as to quantify the deviation of the data with the prediction in some proper way. We will accept the hypothesis if the deviation is not abnormally high. We will treat this problem in the chapter goodness-of-fit tests.

In our third example the number of hypotheses is infinite. To each value of the unknown parameter, i.e. to each di erent mean life, corresponds a di erent prediction. The di culties are very similar to those in case one. If we want to quote probabilities, we are forced to introduce a priori probabilities – here for the parameter under investigation. Again, in most cases no reliable prior information will be available but we can avoid the subjective part by documenting the results in a sensible way. We

134 6 Parameter Inference I

Fig. 6.1. Quantitative Venn diagram. The areas indicate the probabilities for certain combinations of hypotheses Hi and discrete events of type kj .

We assume that the probability P {k|Hi} to observe k is given for a finite number

Here P {Hi} is the assumed probability for the validity of hypothesis i before the observation happens, it is the a priori probability, in short the prior.

In Fig. 6.1 we illustrate relation (6.2) in form of a so called Venn diagram where in the present example 3 out of the 5 hypotheses have the same prior. Each hypothesis bin is divided into 3 regions with areas proportional to the probabilities to observe k = k1, k = k2 and k = k3, respectively. For example when the observation is k = k2 (shadowed in gray) then the gray areas provide the relative probabilities of the validity of the corresponding hypotheses. In our example hypothesis H3 is the most probable, H1 the most unlikely.

The computation of P {k} which is the marginal distribution of k, i.e. the probability of a certain observation, summed over all hypotheses, yields:

X

P {k} = P {k|Hi}P {Hi} .

i

As required, P {Hi|k} is normalized in such a way that the probability that any of the hypotheses is fulfilled is equal to one. We get

In words: The probability for the validity of hypothesis Hi after the measurement k is equal to the prior P {Hi} of Hi multiplied with the probability to observe k if Hi applies and divided by a normalization factor. When we are only interested in the

Fig. 6.2. Probability density for the true decay time. The mean decay time is 1, the observed value is 1.5.

Example 71. Time of a decay with exponential prior

A detector with finite resolution registers at time t the decay of a K meson. The time resolution corresponds to a Gaussian with variance σ2. We are interested in the time θ at which the decay occurred. The mean lifetime τ of kaons is known. The probability density for the parameter θ before the measurement, the prior, is π(θ) = e−θ/τ /τ, θ ≥ 0. The probability density for t with θ fixed is the Gaussian. Applying (6.5) we obtain the probability density f(θ) = f(θ|t) of the parameter θ,

e−(t−θ)2/(2σ2)e−θ/τ

f(θ) = R ∞ e−(t−θ)2/(2σ2)e−θ/τ dθ ,

0

which is displayed in Fig. 6.2. As a consequence of the exponential prior it is visibly shifted to the left with respect to the observation.

If the value of the probability density fx(x|θ) in (6.5) varies much more rapidly with θ than the prior – this is the case when the observation restricts the parameter drastically – then to a good approximation the prior can be regarded as constant in the interesting region. We then have

fθ(θ|x) ≈ R ∞ . −∞ fx(x|θ)dθ

In this approximation the probability density fθ of the parameter corresponds to the normalized likelihood function.

In practice, fθ often follows to a good approximation a normal distribution. The

decreased by the factor e1/2 define a standard deviation error interval and thus fix

the uncertainty of the estimate ˆ.

θ