double
CalcHalfNormChiLik(const double x, const double n, const double gamma, const double phi)
{
    if (x < DBL_MIN)
    {
        return(0.0);
    }
    else
    {
         double logp = (n-1.0) * log(x) - (0.5 * phi * x * x) + (gamma * x);
         return(exp(logp));
//        return(pow(x, n-1.0) * exp((-0.5 * phi * x * x) + (gamma * x)));
    }
}


/* Calculates the normalizing constant for the scale factor PDF:

   P(x) \propto  x^(n-1) e^-(phi/2 x^2 - gamma x)

   The integral for this can be found in Gradshteyn and Ryzhik,
   p. 365, formula 3.462(1).
*/
double
CalcNormConst(const double n, const double gamma, const double phi)
{
    double      tmpx;

     tmpx = (pow(phi, -0.5 * n) * exp(gamma*gamma / (4.0 * phi))) *
            (tgamma(n) * CalcDnz(-n, -gamma / sqrt(phi)));

//        tmpx = (pow(phi, -0.5 * n) * exp(gamma*gamma / (4.0 * phi))) *
//           (tgamma(n) * CalcUab_large_a(n-0.5, -gamma / sqrt(phi)));

    return(1.0/tmpx);
}


static double
CalcNormConstMm(const double n, const double gamma, const double phi)
{
    double      tmpx;

    tmpx = pow(2.0, 0.5*(n-3.0)) * pow(phi,-0.5*(n+1.0))
         * (
            sqrt(2.0*phi) * tgamma(0.5*n) * gsl_sf_hyperg_1F1(0.5*n, 0.5, 0.5*gamma*gamma/phi)
            + 2.0 * gamma * tgamma(0.5*(n+1.0)) * gsl_sf_hyperg_1F1 (0.5*(n+1.0), 1.5, 0.5*gamma*gamma/phi)
            );

    return(1.0/tmpx);
}


double
CalcHalfNormChi(const double x, const double n, const double gamma, const double phi)
{
    return(CalcHalfNormChiLik(x, n, gamma, phi) * CalcNormConstMm(n, gamma, phi));
}


double
ExpectScale(const double n, const double gamma, const double phi)
{
    return((n+1.0) * CalcUax(n+0.5, -gamma/sqrt(phi))/(sqrt(phi)*CalcUax(n-0.5,-gamma/sqrt(phi))));
    //return((n+1.0) * CalcDnz(-n-1.0, -gamma/sqrt(phi))/(sqrt(phi)*CalcDnz(-n,-gamma/sqrt(phi))));
}


    /* Newton and Raftery 1994
       Approximate Bayesian inference with the weighted likelihood bootstrap (with discussion).
       Journal of the Royal Statistical Society, Series B, 56:3-48.
       Equation 16, p 22

       They suggest delta = 0.01, something "small".
       Curiously, delta = 0.5 results in the average posterior log likelihood
       (which is different from the average posterior likelihood)
 */
    int cnt = 0;
    double term, fac, delta, num, denom, oldbf, bf = 0.0;

    delta = 0.5;
    term = delta * (burn-burnin) / (1.0 - delta);
    do
    {
        ++cnt;
        oldbf = bf;
        num = denom = 0.0;
        for (i = burnin; i < burn; ++i)
        {
            liksi = liks[i];
            diff = liksi - blik;
            ediff = exp(diff);
            fac = delta*bf + (1.0-delta)*ediff;
            num += ediff / fac;
            denom += 1.0 / fac;
        }

        bf = (term + num) / (term * bf + denom);
        //printf("Marginal likelihood2[%3d]: % 14.2f \n", cnt, log(bf) + blik);
    }
    while (fabs(oldbf - bf) > bf * 1e-7 && cnt < 1000);
    printf("    Marginal likelihood2[%3d]: % 14.6f \n\n", cnt, log(bf) + blik);

////////////////////////////////////////////////////////////////////////////////////////////////////

double
f(double x, void *params)
{
    double *p = (double *) params;
    double n = p[0];
    double gamma = p[1];
    double phi = p[2];
    double r = p[3];
    double f = pow(x,r) * CalcHalfNormChiLik(x, n, gamma, phi);
    //printf("x: %e  n: %e  gamma: %e  phi: %e  prob: %e\n", x, n, gamma, phi, f);
    //fflush(NULL);
    return f;
}



    double Cm, lik;
    double x,n,phi,gamma;
    FILE *metfp = fopen("metropolis.txt", "w");
    FILE *rejfp = fopen("rejection.txt", "w");

    int samples = 100000;

    n=100;
    phi=100;
    gamma=100;

    double scalemax = ScaleMax(n,gamma,phi);
    printf("Mx: % e\n", scalemax);
    printf("FI: % e\n", 1.0 / (phi + (n-1.0)/(scalemax*scalemax)));
    printf("~Ex: % e\n", 1.0+(gamma/sqrt(phi)));
    printf("~Ex_dlt: % e\n", sqrt(n/phi));


    double sm = ScaleMax(n,gamma,phi);
    double width = sqrt(1.0 / (phi + (n-1.0)/(sm*sm)));
    double xv = sm;
    int skip = 3;

    for (i = 0; i < samples; ++i)
    {
        xv = scale_log_met3(n, gamma, phi, xv, r2, sm, width, skip);
        fprintf(metfp, "%e\n", xv);
    }

if (0)
    for (i = 0; i < samples; ++i)
    {
        fprintf(rejfp, "%e\n", scale_rejection(n, gamma, phi, r2));
    }

    fclose(metfp);
    fclose(rejfp);

//     printf("1F1:\n");
//     fflush(NULL);
//     printf("1F1: %e\n", gsl_sf_hyperg_1F1(0.5*(n+1.0), 1.5, 0.5*gamma*gamma/phi));
//     fflush(NULL);

    //C = CalcNormConst(n, gamma, phi);

    if (0)
    {
        Cm = CalcNormConstMm(n, gamma, phi);
        for (i=0;i<100;++i)
        {
            x=i*0.1;
            lik = CalcHalfNormChi(x, n, gamma, phi);

            //printf("L-[%3d]: %f % e\n", i, x, lik);
            //printf("CL[%3d]: %f % e\n", i, x, C*lik);
            printf("CLm[%3d]: %f % e\n", i, x, lik);
        }
    }

    //double integral1 = integrate_romberg_f3(CalcHalfNormChi, n, gamma, phi, 0.0, 100.0);
    //printf("integral: %e\n", integral1);

    //printf("Ex: %e\n", ExpectScale(n,gamma,phi));

    //printf("Dnz[]: % 7.4e\n", CalcDnz(-n, -gamma / sqrt(phi)));
    //printf("Dnz[]: % 7.4e\n", CalcDnz(2, 2));
    //printf("C-: % .18f\n", C);
    Cm = CalcNormConstMm(n, gamma, phi);
    printf("Cm: % .18f\n", Cm);

/*  for (i=0;i<50;++i) */
/*      printf("Uax[-1.5, %f]: % 7.4e\n", i*0.1, CalcUax(5, i*0.1)); */

    //printf("U(): %e\n", gsl_sf_hyperg_U(-2.0, 0.5, 1.13));

////////////////////////////////////////////////////////////////////////////////////////////////////

    #include <gsl/gsl_integration.h>
    gsl_integration_workspace      *w = gsl_integration_workspace_alloc(1000);
    double                         *params = malloc(4 * sizeof(double));
    double                          result, error;
    double                          p, m1, m2, m3, m4, v, s, k, sd;

    params[0] = n;
    params[1] = gamma;
    params[2] = phi;

    gsl_function F;
    F.function = &f;
    F.params = &params[0];

    //gsl_integration_qag(&F, 0.0, 10.0, 0.0, 1e-7, 1000, GSL_INTEG_GAUSS61, w, &result, &error);
    params[3] = 0.0;
    gsl_integration_qagiu(&F, 0.0, 0.0, 1e-6, 1000, w, &result, &error);
    printf ("result (C)      = % .18e +/- % .18e\n", result, error);
    printf ("result (1/C)    = % .18e +/- % .18e\n", 1.0/result, error/(result*result));
    p = 1.0/result;

    params[3] = 1.0;
    gsl_integration_qagiu(&F, 0.0, 0.0, 1e-6, 1000, w, &result, &error);
    printf ("result (m1)     = % .18f +/- % .18f\n", p*result, p*error);
    m1 = p*result;

    params[3] = 2.0;
    gsl_integration_qagiu(&F, 0.0, 0.0, 1e-6, 1000, w, &result, &error);
    printf ("result (m2)     = % .18f +/- % .18f\n", p*result, p*error);
    m2 = p*result;

    params[3] = 3.0;
    gsl_integration_qagiu(&F, 0.0, 0.0, 1e-6, 1000, w, &result, &error);
    printf ("result (m1)     = % .18f +/- % .18f\n", p*result, p*error);
    m3 = p*result;

    params[3] = 4.0;
    gsl_integration_qagiu(&F, 0.0, 0.0, 1e-6, 1000, w, &result, &error);
    printf ("result (m2)     = % .18f +/- % .18f\n", p*result, p*error);
    m4 = p*result;

    v = m2 - m1*m1;
    sd = sqrt(v);
    s = (2.0*m1*m1*m1 - 3.0*m1*m2 + m3)/(sd*sd*sd);
    k = 3.0-(-3.0*m1*m1*m1*m1 + 6.0*m1*m1*m2 - 4.0*m1*m3 + m4)/(sd*sd*sd*sd);
    printf ("emp exp = % e\n", m1);
    printf ("emp var = % e\n", v);
    printf ("emp skw = % e\n", s);
    printf ("emp kur = % e\n", k);

    free(params);
    gsl_integration_workspace_free(w);

////////////////////////////////////////////////////////////////////////////////////////////////////


static double
gamma_large_dev3(const double a, gsl_rng *r2)
{
    /* Works only if a > 1, and is most efficient if a is large

       This algorithm, reported in Knuth, is attributed to Ahrens.
       A faster one, we are told, can be found in: J. H. Ahrens and
       U. Dieter, Computing 12 (1974) 223-246.  */

    double          sqa, x, y, v;

    sqa = sqrt((2.0 * a) - 1.0);

    do
    {
        do
        {
            y = tan(MY_PI * gsl_rng_uniform(r2));
            x = (sqa * y) + a - 1.0;
        }
        while (x <= 0.0);

        v = gsl_rng_uniform(r2);
    }
    while (v > (1.0 + y*y) * exp((a - 1.0) * log(x / (a - 1.0)) - (sqa * y)));

    return (x);
}


static double
gamma_int_dev3(const unsigned int a, gsl_rng *r2)
{
    unsigned int    i;
    double          prod;

    if (a < 12)
    {
        prod  = 1.0;
        for (i = 0; i < a; ++i)
            prod *= gsl_rng_uniform(r2);

        /* Note: for 12 iterations we are safe against underflow, since
           the smallest positive random number is O(2^-32). This means
           the smallest possible product is 2^(-12*32) = 10^-116 which
           is within the range of const double precision. */

        return (-log(prod));
    }
    else
    {
        return (gamma_large_dev3(a, r2));
    }
}


static double
gamma_frac_dev3(const double a, gsl_rng *r2)
{
    /* This is exercise 16 from Knuth; see page 135, and the solution is
       on page 551.  */
    double          p, q, x, u, v;

    p = MY_E / (a + MY_E);

    do
    {
        u = gsl_rng_uniform(r2);
        v = gsl_rng_uniform(r2);

        if (u < p)
        {
            x = exp((1.0 / a) * log(v));
            q = exp(-x);
        }
        else
        {
            x = 1.0 - log(v);
            q = exp((a - 1.0) * log(x));
        }
    }
    while (gsl_rng_uniform(r2) >= q);

    return (x);
}


static double
gamma_dev3(const double b, const double c, gsl_rng *r2)
{
   /* assume a > 0 */
   int nc = floor (c);

   if (c == nc)
       return (b * gamma_int_dev3(nc, r2));
   else if (nc == 0.0)
       return (b * gamma_frac_dev3(c, r2));
   else
       return (b * (gamma_int_dev3(nc, r2) + gamma_frac_dev3(c - nc, r2)));
}


static double
invgamma_dev3(const double b, const double c, gsl_rng *r2)
{
    return(1.0 / gamma_dev3(1.0/b, c, r2));
}