As you are "only" integrating a one dimensional function, have you also tried the standard R function integrate? 

It is well known that numerical integration can be tricky if the integrand does not "behave well". If there is a problem with the integrand, function integrate might give a warning.

You could also check your integrand graphically for different values of your parameters. Is the function fast changing at certain positions or does it show sharp peaks? Maybe, you have to split up the integration interval.

Best

Matthias

Thanks Matthias,

I have used the standard function "integrate". In fact, it gives the same result as function integral when it works (since it fails sometimes whereas "integral" does not). 

I also have ploted the integrand. It decreases fast and does not have picks (it is smooth). It starts with value at 0 equals to 0.1 and decreases at 10000 to a value of 10^(-7) and continues to decrease. I know that the important area is between 0 and 1000 and the remaining area is very small. Still, the problem resides in estimating the remaining area till infinity. Function int gives a value of 0.95 for the whole line (0,Inf) whereas "integral" gives a value of 0.99. Unfortunately for me is that this small difference creates a great difference in the optimum of the objective function.

I have also tried out to split the function into two parts. I'm not sure which function to use in the bounded interval and which one to use for the infinit part. 

In a different experiment on a different objective function, a similar problem occured. Function integrate failed to calculate the integral between 0 and Infinity. I have used integrate between 0 and 10000 and then used function int between 10000 and infinity. 

I'm concerned about whether you may have encountered a similar problem and was in a position to choose between different integration functions.

Diaa.

Dear Diaa,

what is the error/warning message of integrate?

I had also problems with integrate which was the reason why I added function distrExIntegrate to my package distrEx. The function first tries integrate, if integrate fails it switches to function GLIntegrate performing a Gauss-Legendre quadrature over a finite interval.

My problem was different as integration was part of an iterative algorithm and slightly wrong results were no problem. I just wanted to avoid that the integrate error leads to an error (and hence abort) of my iterative algorithm.

Have you investigated where the difference between int and integral comes from? How are the values for both functions for the intervals (0, 1000) and (1000, Inf)?

Best

Matthias

Indeed, your function distrExIntegrate is better than integrate. I'm also using it. There are still cases when it fails to treat infinity cases perfectly and gives an error that it is divergent although it should not. 

I'm also working with an iterative algorithm, and a priori small errors should not count, because they would vanish in the next iterations. However, I feel that the small differences in the calculus of the integral between "distrExIntegrate" and "int" changes highly the outcome. More precisely, with distrExIntegrate, the optimum is at phi_0. With function "int" it is at phi_1. And if I calculate the objective function in the first case, I really get an infimum at phi_0 and the value at phi_1 is greater. And for the second case the same scenario happens.

I have investigated the differences between int and integral (and distrExIntegrate which gives similar results to integral). Functions "int" gives a value of 0.55 whereas functions integral and distrExIntegrate give 0.95. For the whole line,  int gives 0.95 whereas distrExIntegrate gives 0.99. I have a feeling that function "int" underestimate the value of the integral. You'll find in the attatchement a plot of the objective function between 0 and 100.

It's seems to be a really tricky situation. Some more ideas:

1) Are there special values of the parameter where you could solve the integral analytically (maybe using some computer algebra system) and compare these results with the numerical results?

2) Maybe you can identify dominating terms in your integrand and use this to find an optimal split for the interval.

3) Based on 2) it might also be possible to get a hand on the integrand for x to infty.

Would it be possible for you to provide the crucial code (mainly integrand) to get an even better idea of the problem?

Best

Matthias

In fact, I've already tried out Mathematica. I've simplified the integral with some new parametrization and simple change of variables. Mathematica gave me a formula using the gamma and the incomplete gamma function (which is still an integral function). Supposing that Mathematica can evaluate these functions with a great precision, I should prefere the result of function integral (and your function). 

My integrand is merly the squared root of the multimplication of two GPDs with parameters (localtion=0,scale1,shape1) and (localtion=0,scale2,shape2).

My objective function is the difference between two terms basiclly. An integral and a sum. I have a double optimization. The interior is held upon the first set of parameters (maximization), then the exterior one (minimization of the resulting function) over the second set of parameters.

# The objective function

DualFun = function(alpha,phi)

{

part1 = int(function(x) (1e-100+(1+phi[1]*x/phi[2])^(-1-1/phi[1])/phi[2])^gamma*(1e-100+(1+alpha[1]*x/alpha[2])^(-1-1/alpha[1])/alpha[2])^(-gamma+1), Binf, Inf)

part2 = mean(( (1e-100+(1+phi[1]*Y/phi[2])^(-1-1/phi[1])/phi[2])/(1e-100+(1+alpha[1]*Y/alpha[2])^(-1-1/alpha[1])/alpha[2]) )^gamma)

res = (part1)/(gamma-1) - part2/gamma - 1/(gamma*(gamma-1))

res

}

# Interior optimization

Dphi = function(phi)

{

alphaOpt = constrOptim(theta = phi, f = function(alpha){DualFun(alpha,phi)},ui = ui, ci = ci, method = "Nelder-Mead",control = list(fnscale = -1, maxit = 100,reltol=1e-3))

DualFun(alphaOpt$par,phi)

}

# Exterior optimization

phiNew = constrOptim(theta = phiInit, f = Dphi,ui = ui, ci = ci, method = "Nelder-Mead", control = list(maxit = 100, reltol = 1e-3))

My concern is that the same phenomena repeats itself for a 100 experiments. This optimization calculates the so called (MDphiDE) Minimum dual phi divergence estimator (for the Hellinger case). I'm studying the robustness properties so that vector Y contains a contaminated GPD (location=0, scale=3,shape=0.7) sample. The MDphiDE is considered as an approximation of phi-divergences and should be robust. The resulting estimator when using function "int" is robust whereas the result using distrExIntegrate is not robust and gives the same thing as the MLE.

Thanks alot for your discussion.

Diaa.

Some comments/questions:

- Are gamma and Binf additional fixed parameters? Better don't use global variables but pass them through as additional parameters of the functions.

- You probably don't need the 1/(gamma*(gamma-1)) in DualFun but it shouldn't be a problem

- We are currently developing a package RobExtremes which will also include robust estimators for GPD and GEV; see https://r-forge.r-project.org/projects/robast/

If you are not familiar with SVN, I could also send you development versions of the respective packages (better offline).

- We use a quantile trick for computing integrals with GPD densities which restricts the integration to [0,1]. Maybe this could be of help for you too. The respective implementation is at https://r-forge.r-project.org/scm/viewvc.php/branches/robast-1.0/pkg/RobExtremes/R/Expectation.R?view=markup&root=robast

Best

Matthias

I'm sorry for not giving these constants. gamma = 0.5 fixed. Binf = 0 also fixed.

Ofcourse, the term 1/(gamma*(gamma-1)) is not necessary, but I'm keeping it to give the objective function a positive value and hence as it approaches zero, I'm near the optimum.

I'm aware of the existence of many robust estimators, but the subject of my PhD is about this kind of estimators (although I'm begining to believe that they are not robust !).  Thanks for the extra info.

In fact, I do not know what SVN is and it does interest me to know at least what is it.

Finally, I could not understand your integration method. Are you estimating a quantile after which the pdf becomes below a given precision, and then rescaling the whole thing to [0,1] with any classical method?

Diaa.

SVN is short for subversion a version control system; see http://en.wikipedia.org/wiki/Apache_Subversion.

The quantile trick works as follows:

Instead of integrating the integrand with respect to the GPD density, you can integrate the integrand evaluated at the quantiles of the GPD using the density of the uniform distribution; e.g. in R code

f1

Thanks for the extra info. 

In fact, the quantile trick is useful for my case. I've just used it with function "int" and got the same result as with function integrate (with and without the trick). This means that function "int" does not work properly. 

I've never though about using such a trick although very simple, and in fact all my integrals are of the form \int_{R} (p(x)/q(x))^(gamma-1) * p(x)dx. This means that I can always use the change of variable x= P^(-1) (y)  as along as P (the cdf of p) is strictly increasing (one to one fun) and change the integration into a an integral over (0,1).

I'm very grateful for the fruitful discussion and for your great help.

Thanks Matthias.

Diaa.