Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
817 lines (707 sloc) 22.6 KB
/* Statistical routines for mixtures of generalized extreme value distributions.
*
* Contents:
* 1. The ESL_MIXGEV object
* 2. Evaluating densities and distributions
* 3. Generic API routines: for general interface w/ histogram module
* 4. Dumping plots to xmgrace XY format
* 5. Sampling
* 6. ML fitting to complete data
* 7. Example
*
* Xrefs:
* SRE:STL9/139 original implementation
*
* To-do:
* - Fit*() functions should return eslEINVAL on n=0, eslENORESULT
* on failure due to small n. Compare esl_gumbel. xref J12/93.
* SRE, Wed Nov 27 11:02:14 2013
*/
#include "esl_config.h"
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "easel.h"
#include "esl_dirichlet.h" /* for uniform sampling of a probability vector */
#include "esl_gev.h"
#include "esl_minimizer.h"
#include "esl_random.h"
#include "esl_stats.h"
#include "esl_vectorops.h"
#include "esl_mixgev.h"
/****************************************************************************
* 1. The ESL_MIXGEV object
****************************************************************************/
/* Function: esl_mixgev_Create()
*
* Purpose: Creates an object to hold parameters for a <K>-component
* mixture of generalized extreme value distributions.
*
* Parameters in the object are initialized ($q_k =
* \frac{1}{K}$, $\mu_k = 0$, $\lambda_k = 1$, $\alpha_k =
* 0$), but the caller will want to set these according to
* its own purposes.
*
* After an object is created, the caller can constrain any
* of the components to be a Gumbel (that is, constrain
* $\alpha_k = 0$ by calling <esl_mixgev_ForceGumbel(obj,
* k)>.
*
* Args: K - number of components in the mixture
*
* Returns: ptr to newly allocated/initialized <ESL_MIXGEV> object.
*
* Throws: <NULL> on allocation failure.
*/
ESL_MIXGEV *
esl_mixgev_Create(int K)
{
ESL_MIXGEV *mg = NULL;
int k;
int status;
ESL_ALLOC(mg, sizeof(ESL_MIXGEV));
mg->q = mg->mu = mg->lambda = mg->alpha = mg->wrk = NULL;
mg->isgumbel = NULL;
mg->K = K;
ESL_ALLOC(mg->q, sizeof(double) * K);
ESL_ALLOC(mg->mu, sizeof(double) * K);
ESL_ALLOC(mg->lambda, sizeof(double) * K);
ESL_ALLOC(mg->alpha, sizeof(double) * K);
ESL_ALLOC(mg->wrk, sizeof(double) * K);
ESL_ALLOC(mg->isgumbel, sizeof(int) * K);
for (k = 0; k < K; k++)
{
mg->q[k] = 1. / (double) K;
mg->mu[k] = 0.;
mg->lambda[k] = 1.;
mg->alpha[k] = 0.;
mg->isgumbel[k] = FALSE;
}
return mg;
ERROR:
esl_mixgev_Destroy(mg);
return NULL;
}
/* Function: esl_mixgev_Destroy()
*
* Purpose: Deallocates the mixture GEV parameter object <mg>.
*
* Args: mg - ptr to the object to be deallocated.
*
* Returns: (void)
*/
void
esl_mixgev_Destroy(ESL_MIXGEV *mg)
{
if (mg == NULL) return;
if (mg->q != NULL) free(mg->q);
if (mg->mu != NULL) free(mg->mu);
if (mg->lambda != NULL) free(mg->lambda);
if (mg->alpha != NULL) free(mg->alpha);
if (mg->wrk != NULL) free(mg->wrk);
if (mg->isgumbel != NULL) free(mg->isgumbel);
free(mg);
}
/* Function: esl_mixgev_Copy()
*
* Purpose: Makes a copy of the mixture GEV parameter object <src>
* in <dest>. Caller must have already allocated <dest> to have
* (at least) the same number of components as <src>.
*
* Args: src - object to be copied
* dest - allocated object to copy <src> into
*
* Returns: <eslOK> on success.
*
* Throws: <eslEINCOMPAT> if <dest> isn't allocated with enough
* components to hold a copy of <src>.
*/
int
esl_mixgev_Copy(ESL_MIXGEV *src, ESL_MIXGEV *dest)
{
int k;
if (dest->K < src->K)
ESL_EXCEPTION(eslEINCOMPAT, "mixture GEV too small to copy into");
for (k = 0; k < src->K; k++)
{
dest->q[k] = src->q[k];
dest->mu[k] = src->mu[k];
dest->lambda[k] = src->lambda[k];
dest->alpha[k] = src->alpha[k];
dest->isgumbel[k] = src->isgumbel[k];
}
dest->K = src->K;
return eslOK;
}
/* Function: esl_mixgev_ForceGumbel()
*
* Purpose: Constrain component <which> of the mixture GEV <mg>
* to be a Gumbel (that is, constrain $\alpha=0$ for
* that component. This constraint will be obeyed by
* any subsequent calls to parameter fitting routines.
*
* Normally would be called just after creating the <mg>
* object, as part of its configuration before trying to
* fit some observed data to a mixture GEV.
*
* Args: mg - mixture GEV object being configured
* which - which component to constrain to a Gumbel
*
* Returns: <eslOK> on success.
*/
int
esl_mixgev_ForceGumbel(ESL_MIXGEV *mg, int which)
{
mg->isgumbel[which] = TRUE;
return eslOK;
}
/*----------------- end ESL_MIXGEV object maintenance ----------------------*/
/****************************************************************************
* 2. Evaluating densities and distributions
****************************************************************************/
/* Function: esl_mixgev_pdf()
*
* Purpose: Returns the probability density function $P(X=x)$ for
* quantile <x>, given mixture GEV parameters <mg>.
*/
double
esl_mixgev_pdf(double x, ESL_MIXGEV *mg)
{
double pdf = 0.;
int k;
for (k = 0; k < mg->K; k++)
pdf += mg->q[k] * esl_gev_pdf(x, mg->mu[k], mg->lambda[k], mg->alpha[k]);
return pdf;
}
/* Function: esl_mixgev_logpdf()
*
* Purpose: Returns the log of the PDF ($\log P(X=x)$) for quantile <x>,
* given mixture GEV parameters <mg>.
*/
double
esl_mixgev_logpdf(double x, ESL_MIXGEV *mg)
{
int k;
for (k = 0; k < mg->K; k++)
if (mg->q[k] == 0.0)
mg->wrk[k] = -eslINFINITY;
else
mg->wrk[k] = log(mg->q[k]) +
esl_gev_logpdf(x, mg->mu[k], mg->lambda[k], mg->alpha[k]);
return esl_vec_DLogSum(mg->wrk, mg->K);
}
/* Function: esl_mixgev_cdf()
*
* Purpose: Returns the cumulative distribution function $P(X \leq x)$
* for quantile <x>, given mixture GEV parameters <mg>.
*/
double
esl_mixgev_cdf(double x, ESL_MIXGEV *mg)
{
double cdf = 0.;
int k;
for (k = 0; k < mg->K; k++)
cdf += mg->q[k] * esl_gev_cdf(x, mg->mu[k], mg->lambda[k], mg->alpha[k]);
return cdf;
}
/* Function: esl_mixgev_logcdf()
*
* Purpose: Returns the log of the CDF $\log P(X \leq x)$
* for quantile <x>, given mixture GEV parameters <mg>.
*/
double
esl_mixgev_logcdf(double x, ESL_MIXGEV *mg)
{
int k;
for (k = 0; k < mg->K; k++)
if (mg->q[k] == 0.0)
mg->wrk[k] = -eslINFINITY;
else
mg->wrk[k] = log(mg->q[k]) +
esl_gev_logcdf(x, mg->mu[k], mg->lambda[k], mg->alpha[k]);
return esl_vec_DLogSum(mg->wrk, mg->K);
}
/* Function: esl_mixgev_surv()
*
* Purpose: Returns the survivor function $P(X > x)$ (1-CDF)
* for quantile <x>, given mixture GEV parameters <mg>.
*/
double
esl_mixgev_surv(double x, ESL_MIXGEV *mg)
{
double srv = 0.;
int k;
for (k = 0; k < mg->K; k++)
srv += mg->q[k] * esl_gev_surv(x, mg->mu[k], mg->lambda[k], mg->alpha[k]);
return srv;
}
/* Function: esl_mixgev_logsurv()
*
* Purpose: Returns the log survivor function $\log P(X > x)$ (log(1-CDF))
* for quantile <x>, given mixture GEV parameters <mg>.
*/
double
esl_mixgev_logsurv(double x, ESL_MIXGEV *mg)
{
int k;
for (k = 0; k < mg->K; k++)
{
mg->wrk[k] = log(mg->q[k]);
mg->wrk[k] += esl_gev_logsurv(x, mg->mu[k], mg->lambda[k], mg->alpha[k]);
}
return esl_vec_DLogSum(mg->wrk, mg->K);
}
/* Function: esl_mixgev_invcdf()
*
* Purpose: Calculates the inverse CDF for a mixture GEV <mg>,
* returning the quantile <x> at which the CDF is <p>,
* where $0 < p < 1$.
*
* The inverse CDF of a mixture model has no analytical
* expression as far as I'm aware. The calculation here is
* a brute force bisection search in <x> using the CDF
* function. It will suffice for a small number of calls
* (for plotting applications, for example), but beware, it is not
* efficient.
*/
double
esl_mixgev_invcdf(double p, ESL_MIXGEV *mg)
{
double x1, x2, xm; /* low, high guesses at x */
double f1, f2, fm;
double tol = 1e-6;
x2 = esl_vec_DMin(mg->mu, mg->K);
x1 = x2 - 1.;
do { /* bracket, left side */
x1 = x1 + 2.*(x2-x1);
f1 = esl_mixgev_cdf(x1, mg);
} while (f1 > p);
do { /* bracket, right side */
x2 = x2 + 2.*(x2-x1);
f2 = esl_mixgev_cdf(x2, mg);
} while (f2 < p);
do { /* bisection */
xm = (x1+x2) / 2.;
fm = esl_mixgev_cdf(xm, mg);
if (fm > p) x2 = xm;
else if (fm < p) x1 = xm;
else return xm; /* unlikely case of fm==p */
} while ( (x2-x1)/(x1+x2+1e-9) > tol);
xm = (x1+x2) / 2.;
return xm;
}
/*-------------------- end densities & distributions ------------------------*/
/****************************************************************************
* 3. Generic API routines: for general interface w/ histogram module
****************************************************************************/
/* Function: esl_mixgev_generic_pdf()
*
* Purpose: Generic-API wrapper around <esl_mixgev_pdf()>, taking
* a void ptr to a <ESL_MIXGEV> parameter structure.
*/
double
esl_mixgev_generic_pdf(double x, void *params)
{
ESL_MIXGEV *mg = (ESL_MIXGEV *) params;
return esl_mixgev_pdf(x, mg);
}
/* Function: esl_mixgev_generic_cdf()
*
* Purpose: Generic-API wrapper around <esl_mixgev_cdf()>, taking
* a void ptr to a <ESL_MIXGEV> parameter structure.
*/
double
esl_mixgev_generic_cdf(double x, void *params)
{
ESL_MIXGEV *mg = (ESL_MIXGEV *) params;
return esl_mixgev_cdf(x, mg);
}
/* Function: esl_mixgev_generic_surv()
*
* Purpose: Generic-API wrapper around <esl_mixgev_surv()>, taking
* a void ptr to a <ESL_MIXGEV> parameter structure.
*/
double
esl_mixgev_generic_surv(double x, void *params)
{
ESL_MIXGEV *mg = (ESL_MIXGEV *) params;
return esl_mixgev_surv(x, mg);
}
/* Function: esl_mixgev_generic_invcdf()
*
* Purpose: Generic-API wrapper around <esl_mixgev_invcdf()>, taking
* a void ptr to a <ESL_MIXGEV> parameter structure.
*/
double
esl_mixgev_generic_invcdf(double p, void *params)
{
ESL_MIXGEV *mg = (ESL_MIXGEV *) params;
return esl_mixgev_invcdf(p, mg);
}
/*------------------------ end generic API ---------------------------------*/
/****************************************************************************
* 4. Dumping plots to xmgrace XY format
****************************************************************************/
/* Function: esl_mixgev_Plot()
*
* Purpose: Plot some function <func> (for instance, <esl_mixgev_pdf()>)
* for mixture GEV parameters <mg>, for a range of
* quantiles x from <xmin> to <xmax> in steps of <xstep>;
* output to an open stream <fp> in xmgrace XY input format.
*
* Returns: <eslOK> on success.
*
* Throws: <eslEWRITE> on any system write error.
*/
int
esl_mixgev_Plot(FILE *fp, ESL_MIXGEV *mg,
double (*func)(double x, ESL_MIXGEV *mg),
double xmin, double xmax, double xstep)
{
double x;
for (x = xmin; x <= xmax; x += xstep)
if (fprintf(fp, "%f\t%g\n", x, (*func)(x, mg)) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "mixgev plot write failed");
if (fprintf(fp, "&\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "mixgev plot write failed");
return eslOK;
}
/*-------------------- end plot dumping routines ---------------------------*/
/****************************************************************************
* 5. Sampling
****************************************************************************/
/* Function: esl_mixgev_Sample()
*
* Purpose: Sample a random variate x from a mixture GEV <mg>,
* given random number source <r>.
*/
double
esl_mixgev_Sample(ESL_RANDOMNESS *r, ESL_MIXGEV *mg)
{
int k;
k = esl_rnd_DChoose(r, mg->q, mg->K);
return esl_gev_Sample(r, mg->mu[k], mg->lambda[k], mg->alpha[k]);
}
/*--------------------------- end sampling ---------------------------------*/
/****************************************************************************
* 6. ML fitting to complete data
****************************************************************************/
struct mixgev_data {
double *x;
int n;
double *wrk; /* workspace vector [0..K-1]*/
ESL_MIXGEV *mg;
};
/* Given mixture GEV parameters in <mg>;
* do appropriate c.o.v.'s to unconstrained real parameters
* and fill in the packed parameter vector <p>.
*
* First K-1 are Q_1..Q_K-1 mixture coefficient parameters; Q_0 implicitly 0;
* cov is q_k = e^{Q_k} / \sum_j e^{Q_j},
* so Q_k = log(q_k) - log(q_0).
* Then K components: mu, lambda, optional alpha;
* mu, alpha are already unconstrained real;
* lambda cov is lambda = e^w, w = log(lambda).
*/
static void
mixgev_pack_paramvector(double *p, int np, ESL_MIXGEV *mg)
{
int i; /* counter in parameter vector p */
int k; /* counter in mixture components */
double z; /* tmp variable */
/* mixture coefficients */
z = log(mg->q[0]);
i = 0;
for (k = 1; k < mg->K; k++)
p[i++] = log(mg->q[k]) - z;
/* gev parameters */
for (k = 0; k < mg->K; k++)
{
p[i++] = mg->mu[k];
p[i++] = log(mg->lambda[k]);
if (! mg->isgumbel[k]) p[i++] = mg->alpha[k];
}
/* assert(i==np) in debugging, if you want */
}
/* Same as above but in reverse: given parameter vector <p>,
* do appropriate c.o.v. back to desired parameter space, and
* fill in the mixture GEV structure <mg>.
*/
static void
mixgev_unpack_paramvector(double *p, int np, ESL_MIXGEV *mg)
{
int i; /* counter in parameter vector p */
int k; /* counter in mixture components */
double z; /* tmp variable */
/* Fetch the params in their c.o.v. space first
*/
i = 0;
mg->q[0] = 0; /* implicitly */
for (k = 1; k < mg->K; k++)
mg->q[k] = p[i++];
for (k = 0; k < mg->K; k++)
{
mg->mu[k] = p[i++];
mg->lambda[k] = p[i++];
if (!mg->isgumbel[k]) mg->alpha[k] = p[i++];
else mg->alpha[k] = 0.;
}
assert(i==np);
/* Convert mix coefficients back to probabilities;
* their c.o.v. is q_k = e^{Q_k} / \sum_k e^{Q_k}
* which rearranges to exp(Q_k - log[\sum_k e^Q_k]),
* and we have the DLogSum() function to compute the log sum.
*/
z = esl_vec_DLogSum(mg->q, mg->K);
for (k = 0; k < mg->K; k++)
mg->q[k] = exp(mg->q[k] - z);
/* lambda c.o.v. is \lambda = e^w
*/
for (k = 0; k < mg->K; k++)
mg->lambda[k] = exp(mg->lambda[k]);
}
static double
mixgev_complete_func(double *p, int np, void *dptr)
{
struct mixgev_data *data = (struct mixgev_data *) dptr;
ESL_MIXGEV *mg = data->mg;
int i;
double logL;
/* Use the current parameter vector (in its unconstrained
* real c.o.v. space) to deduce what the current mixture GEV
* parameters are:
*/
mixgev_unpack_paramvector(p, np, mg);
/* Calculate the log likelihood:
*/
logL = 0;
for (i = 0; i < data->n; i++)
logL += esl_mixgev_logpdf(data->x[i], mg);
/* return the NLL
*/
return -logL;
}
/* Function: esl_mixgev_FitGuess()
*
* Purpose: Make initial randomized guesses at the parameters
* of mixture GEV <mg>, using random number generator
* <r> and observed data consisting of <n> values
* <x[0..n-1]>. This guess is a suitable starting
* point for a parameter optimization routine, such
* as <esl_mixgev_FitComplete()>.
*
* Specifically, we estimate one 'central' guess
* for a single Gumbel fit to the data, using the
* method of moments. Then we add $\pm 10\%$ to that 'central'
* $\mu$ and $\lambda$ to get each component
* $\mu_i$ and $\lambda_i$. The $\alpha_i$ parameters
* are generated by sampling uniformly from $-0.1..0.1$.
* Mixture coefficients $q_i$ are sampled uniformly.
*
* Args: r - randomness source
* x - vector of observed data values to fit, 0..n-1
* n - number of values in <x>
* mg - mixture GEV to put guessed params into
*
* Returns: <eslOK> on success.
*/
int
esl_mixgev_FitGuess(ESL_RANDOMNESS *r, double *x, int n, ESL_MIXGEV *mg)
{
double mean, variance;
double mu, lambda;
int k;
esl_stats_DMean(x, n, &mean, &variance);
lambda = eslCONST_PI / sqrt(6.*variance);
mu = mean - 0.57722/lambda;
esl_dirichlet_DSampleUniform(r, mg->K, mg->q);
for (k = 0; k < mg->K; k++)
{
mg->mu[k] = mu + 0.2 * mu * (esl_random(r) - 0.5);
mg->lambda[k] = lambda + 0.2 * lambda * (esl_random(r) - 0.5);
if (mg->isgumbel[k]) mg->alpha[k] = 0.;
else mg->alpha[k] = 0.2 * (esl_random(r) - 0.5);
}
return eslOK;
}
/* Function: esl_mixgev_FitComplete()
*
* Purpose: Given <n> observed data values <x[0..n-1]>, and
* an initial guess at a mixture GEV fit to those data
* <mg>, use conjugate gradient descent to perform
* a locally optimal maximum likelihood mixture
* GEV parameter fit to the data.
*
* To obtain a reasonable initial guess for <mg>,
* see <esl_mixgev_FitGuess()>.
*
* Args: x - observed data, <x[0..n-1]>.
* n - number of samples in <x>
* mg - mixture GEV to estimate, w/ params set to
* an initial guess.
*
* Returns: <eslOK> on success, and <mg> contains local
* ML estimate for mixture GEV parameters.
*
* Throws: <eslEMEM> on allocation error, and <mg> is unchanged
* from its initial state.
*/
int
esl_mixgev_FitComplete(double *x, int n, ESL_MIXGEV *mg)
{
struct mixgev_data data;
int status;
double *p = NULL;
double *u = NULL;
double *wrk = NULL;
double tol;
int np;
double fx;
int k;
int i;
tol = 1e-6;
/* Determine number of free parameters and allocate
*/
np = mg->K-1; /* K-1 mix coefficients free */
for (k = 0; k < mg->K; k++)
np += (mg->isgumbel[k])? 2 : 3;
ESL_ALLOC(p, sizeof(double) * np);
ESL_ALLOC(u, sizeof(double) * np);
ESL_ALLOC(wrk, sizeof(double) * np * 4);
/* Copy shared info into the "data" structure
*/
data.x = x;
data.n = n;
data.wrk = wrk;
data.mg = mg;
/* From mg, create the parameter vector.
*/
mixgev_pack_paramvector(p, np, mg);
/* Define the step size vector u.
*/
i = 0;
for (k = 1; k < mg->K; k++) u[i++] = 1.0;
for (k = 0; k < mg->K; k++)
{
u[i++] = 1.0;
u[i++] = 1.0;
if (! mg->isgumbel[k]) u[i++] = 0.02;
}
ESL_DASSERT1( (np == i) );
/* Feed it all to the mighty optimizer.
*/
status = esl_min_ConjugateGradientDescent(p, u, np, &mixgev_complete_func, NULL,
(void *) (&data), tol, wrk, &fx);
if (status != eslOK) goto ERROR;
/* Convert the final parameter vector back to a mixture GEV
*/
mixgev_unpack_paramvector(p, np, mg);
free(p);
free(u);
free(wrk);
return eslOK;
ERROR:
if (p != NULL) free(p);
if (u != NULL) free(u);
if (wrk != NULL) free(wrk);
return status;
}
/*--------------------------- end fitting ----------------------------------*/
/****************************************************************************
* Example main()
****************************************************************************/
#ifdef eslMIXGEV_EXAMPLE
/*::cexcerpt::mixgev_example::begin::*/
/* compile:
gcc -g -Wall -I. -L. -o example -DeslMIXGEV_EXAMPLE esl_mixgev.c -leasel -lm
* run: ./example
*/
#include <stdio.h>
#include <stdlib.h>
#include "easel.h"
#include "esl_mixgev.h"
#include "esl_random.h"
int
main(int argc, char **argv)
{
FILE *fp;
ESL_RANDOMNESS *r; /* source of random numbers */
ESL_MIXGEV *mg; /* mixture GEV to sample from */
ESL_MIXGEV *emg; /* estimated mixture GEV */
double *x; /* sampled dataset */
int n = 100000; /* number of samples */
int i;
int k;
double nll;
double min, max;
r = esl_randomness_Create(42);
mg = esl_mixgev_Create(2);
mg->q[0] = 0.85; mg->q[1] = 0.15;
mg->mu[0] = -2.72; mg->mu[1] = -2.0;
mg->lambda[0] = 2.5; mg->lambda[1] = 1.0;
mg->alpha[0] = 0.; mg->alpha[1] = 0.09;
nll = 0.;
min = 99999;
max = -99999;
x = malloc(sizeof(double) * n);
for (i = 0; i < n; i++)
{
x[i] = esl_mixgev_Sample(r, mg);
nll -= esl_mixgev_logpdf(x[i], mg);
if (x[i] > max) max = x[i];
if (x[i] < min) min = x[i];
}
printf("NLL of known mixGEV: %g\n", nll);
/* Dump the raw data samples to an R file.
*/
fp = fopen("data.out", "w");
fprintf(fp, " val\n");
for (i = 0; i < n; i++)
fprintf(fp, "%d %f\n", i+1, x[i]);
fclose(fp);
emg = esl_mixgev_Create(2);
esl_mixgev_FitGuess(r, x, n, emg);
/* esl_mixgev_Copy(mg, emg); */
esl_mixgev_ForceGumbel(emg, 0);
esl_mixgev_FitComplete(x, n, emg);
printf("Component q mu lambda alpha\n");
for (k=0; k < 2; k++)
printf("%d\t%7.4f\t%7.2f\t%7.4f\t%7.4f\n",
k, emg->q[k], emg->mu[k], emg->lambda[k], emg->alpha[k]);
nll = 0.;
for (i = 0; i < n; i++)
nll -= esl_mixgev_logpdf(x[i], emg);
printf("NLL of fitted mixGEV: %g\n", nll);
/* Dump some R commands for showing these distributions
*/
printf("library(ismev)\n");
printf("library(evd)\n");
printf("d <- read.table(\"data.out\")$val\n");
printf("plot(density(d,bw=0.2), log=\"y\")\n");
printf("min <- %f\n", min);
printf("max <- %f\n", max);
printf("xax <- seq(min-2, max+5, by=0.1)\n");
printf("cc <- xax - xax\n");
printf("zc <- xax - xax\n");
for (k = 0; k < mg->K; k++)
{
printf("c%d <- %f * dgev(xax, %f, %f, %f)\n",
k, mg->q[k], mg->mu[k], 1./mg->lambda[k], mg->alpha[k]);
printf("cc <- cc + c%d\n", k);
printf("lines(xax, c%d, col=\"blue\")\n", k);
}
for (k = 0; k < emg->K; k++)
{
printf("z%d <- %f * dgev(xax, %f, %f, %f)\n",
k, emg->q[k], emg->mu[k], 1./emg->lambda[k], emg->alpha[k]);
printf("zc <- zc + z%d\n", k);
printf("lines(xax, z%d, col=\"blue\")\n", k);
}
printf("lines(xax, cc, col=\"green\")\n");
printf("lines(xax, zc, col=\"red\")\n");
esl_mixgev_Destroy(mg);
esl_mixgev_Destroy(emg);
esl_randomness_Destroy(r);
free(x);
return 0;
}
/*::cexcerpt::mixgev_example::end::*/
#endif /*eslMIXGEV_EXAMPLE*/
You can’t perform that action at this time.