Download Practical Guide To Computer Simulations

76
1
2
3
4
5
6
A.K. Hartmann: Practical Guide to Computer Simulations
double ll_ft(const gsl_vector *par, void *param)
{
double lambda, x0;
/* parameters of pdf
sample_t *sample;
/* sample
double sum;
/* sum of log-likelihood contributions
int i;
/* loop counter
*/
*/
*/
*/
7
lambda = gsl_vector_get(par, 0);
x0 = gsl_vector_get(par, 1);
sample = (sample_t *) param;
8
9
10
/* get data */
11
sum = sample->n*log(lambda);
/* calculate log likelihood */
for(i=0; i<sample->n; i++)
sum -= lambda*(sample->x[i]-x0) +
exp(-lambda*(sample->x[i]-x0));
12
13
14
15
16
return(-sum);
17
18
/* return - log likelihood */
}
First, we convert the pointers passed as arguments to the data format that
we find useful (lines 8–10). Next, the actual log likelihood
l(λ, xo ) = n log λ − λ
n−1
X
i=0
(xi − x0 ) −
n−1
X
i=0
exp(−λ(xi − x0 ))
is calculated in lines 12–15 and finally returned with inverted sign (line 17).
The GSL has built in several minimization algorithms. They are all put
under one of two frameworks. One framework is for algorithms which require the
target function and its first derivatives. The other framwork contains algorithms
where just the target function is sufficient. Here we use the simplex algorithm,
which belongs to the latter form. It works by spanning a simplex,16 evaluating
the target functions at the corners of the simplex, and iteratively changing the
simplex until it is very small and contains the solution. Note that the algorithm
is only able to find local minima, and only one of them. If several minima exist,
the choice of the initial parameters strongly influence the final results; Here,
one maybe has to try several parameters. For details see [Galassi et al. (2006)].
Here we only show how to use the minimizer. The minimizer itself is stored in a
special data structure of type gsl_multimin_fminimizer. The target function
has to be put into a “surrounding” variable of type gsl_multimin_function.
Furthermore, one needs two gsl_vector variables to store the current estimate
for the optimum (specifying the position of the simplex) and to store the size of
the simplex. Also, par is used here to state the dimension of the target function
argument (2) and sample to store the sample.
These variables are declared as follows:
16 A
simplex is a convex set in an n-dimensional space generated by n + 1 corner points.