Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?

[Martin Jansche]

Most of these look sensible. I'd be inclined to try VEGAS and compare
against plain Monte Carlo.

On Mon, Jan 1, 2018 at 6:10 AM, Vasu Jaganath <address@hidden>
wrote:

> HI Martin, Yes one of my Q is very discontinuous with respect to my
> integration variable Z.
>
> I have attached the plots for you to see
>
> On Sun, Dec 31, 2017 at 9:20 PM, Vasu Jaganath <address@hidden>
> wrote:
>
>> Yes, I will show you the plots soon, Q is actually 2 variable function
>> but for my calculations I am treating one of the variables as a parameter,
>> which is a physically valid assumption. Yes I do encounter alpha and beta
>> values less than 1.
>>
>> On Sun, Dec 31, 2017, 9:13 PM Martin Jansche <address@hidden> wrote:
>>
>>> So you want to find E[f] = \int_0^1 f(x) dbeta(x | a, b) dx. Can you plot
>>> your typical f? And what are typical values of the parameters a and b? Do
>>> you encounter a<=1 or b<=1? If so, how does f(x) behave as x approaches 0
>>> or 1?
>>>
>>> On Mon, Jan 1, 2018 at 3:37 AM, Patrick Alken <address@hidden>
>>> wrote:
>>>
>>> > The question is whether your Q contains any singularities, or is highly
>>> > oscillatory? Is such cases fixed point quadrature generally doesn't do
>>> > well. If Q varies fairly smoothly over your interval, you should give
>>> > fixed point quadrature a try and report back if it works well enough
>>> for
>>> > your problem. The routines you want are documented here:
>>> >
>>> > http://www.gnu.org/software/gsl/doc/html/integration.html#
>>> > fixed-point-quadratures
>>> >
>>> > Also, if QAGS isn't working well for you, try also the CQUAD routines.
>>> > I've found CQUAD is more robust than QAGS in some cases
>>> >
>>> > On 12/31/2017 05:28 PM, Vasu Jaganath wrote:
>>> > > I have attached my entire betaIntegrand function. It is a bit
>>> complicated
>>> > > and very boiler-plate, It's OpenFOAM code (where scalar = double
>>> etc.) I
>>> > > hope you can get the jist from it.
>>> > > I can integrate the Q using monte-carlo sampling integration.
>>> > >
>>> > > Q is nothing but tabulated values of p,rho, mu etc. I lookup Q using
>>> the
>>> > > object "solver" in the snippet.
>>> > >
>>> > > I have verified evaluating <Q> when I am not using those <Q> values
>>> back
>>> > in
>>> > > the solution, It works OK.
>>> > >
>>> > > Please ask me anything if it seems unclear.
>>> > >
>>> > >
>>> > >
>>> > >
>>> > >
>>> > >
>>> > > On Sun, Dec 31, 2017 at 3:55 PM, Martin Jansche <address@hidden>
>>> > wrote:
>>> > >
>>> > >> Can you give a concrete example of a typical function Q?
>>> > >>
>>> > >> On Sat, Dec 30, 2017 at 3:42 AM, Vasu Jaganath <
>>> > address@hidden>
>>> > >> wrote:
>>> > >>
>>> > >>> 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.
>>> > >>>
>>> > >>> Thanks,
>>> > >>> Vasu
>>> > >>>
>>> > >>
>>> >
>>> >
>>> >
>>>
>>
>