[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] fixed point or adaptive integration for calculating moments u

From: Vasu Jaganath
Subject: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?
Date: Fri, 29 Dec 2017 20:42:21 -0700

Hi forum,

I am trying to integrate moments, basically first order moments <Q>, i.e.
averages of some flow fields like temperature, density and mu. I am
assuming they distributed according to beta-PDF which is parameterized on
variable Z, whose mean and variance i am calculating separately and using
it to define the shape of the beta-PDF, I have a cut off of not using the
beta-PDF when my mean Z value, i.e <Z> is less than a threshold.

I am using qags, the adaptive integration routine to calculate the moment
integral, however I am restricted to threshold of <Z> = 1e-2.

It complains that the integral is too slowly convergent. However physically
my threshold should be around 5e-5 atleast.

I can integrate these moments with threshold upto 5e-6, using Monte-Carlo
integration, by generating random numbers which are beta-distributed.

Should I look into fixed point integration routines? What routines would
you suggest?

Here is the (very simplified) code snippet where, I calculate alpha and
beta parameter of the PDF

                    zvar   = min(zvar,0.9999*zvar_lim);
                    alpha = zmean*((zmean*(1 - zmean)/zvar) - 1);
                    beta = (1 - zmean)*alpha/zmean;

                    // inside the fucntion to be integrated
                    // lots of boilerplate for Q(x)
                    f = Q(x) * gsl_ran_beta_pdf(x, alpha, beta);

                   // my integration call

                   helper::gsl_integration_qags (&F, 0, 1, 0, 1e-2, 1000,
                                                  w, &result, &error);

And also, I had to give relative error pretty large, 1e-2. However some of
Qs are pretty big order of 1e6.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]