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

[Patrick Alken]

You might be able to used fixed quadrature - you will probably need a
large number of nodes to capture that sharp feature. I would compare
both the CQUAD and fixed-point methods to see which works better for you.

Patrick

On 12/31/2017 11:10 PM, Vasu Jaganath 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 <mailto: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
>     <mailto: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 <mailto: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#
>         <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 <mailto: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 <mailto: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
>         > >>>
>         > >>
>         >
>         >
>         >
>
>