Understanding and implementing numerical integrati

问题:

I need to calculate this integral below, using R:

The q_theta(x) function I managed to do in R with quantile regression (package: quantreg).

matrix=structure(c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 
0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 
0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 
0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 
0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 
0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 
0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 
0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 
0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 
0.98, 0.99, -22.2830664155772, -22.2830664155772, -19.9298291765612, 
-18.2066426767652, -15.2657135034479, -14.921522915965, -13.5035945028536, 
-13.1557269916064, -12.9495709618481, -11.6168348488161, -11.3999095021713, 
-10.6962766764396, -10.0588239375837, -9.12944363439522, -8.15648778610587, 
-8.04133299299019, -7.66558386420434, -7.50906566627427, -6.95626096568998, 
-6.90630556403136, -6.53374879831376, -6.39324677042686, -6.20705804899049, 
-6.09754765999465, -5.91272058217526, -5.75771166206242, -5.3770131257001, 
-5.20892464393192, -5.07372162687422, -4.96706814289334, -4.64404095131293, 
-4.1567394053577, -4.13209444755342, -3.85483644113723, -3.64855238293205, 
-3.53054113507559, -3.46035383338799, -3.03155417364444, -2.93100183005178, 
-2.90491824855193, -2.64056616049773, -2.51857727614607, -2.25163805172486, 
-2.00934783937474, -1.89925824841417, -1.71405007411747, -1.65905834683964, 
-1.47502511311988, -1.42755073292529, -1.20464216637298, -1.08574103345057, 
-0.701134735371922, -0.590656010656201, -0.290335898959635, -0.0575062007348038, 
0.0778328375033378, 0.165234593185889, 0.230651883848336, 0.316817885358695, 
0.34841775605248, 0.516869604496075, 0.59743162507581, 0.857843937404964, 
0.939734010162078, 1.12533017928147, 1.27037182428776, 1.52040854525927, 
1.76577933448152, 2.07456447851822, 2.17389787235523, 2.27567786362425, 
2.3850323163509, 2.55365596853891, 2.61208242890655, 2.77359226593771, 
2.93275094039929, 3.07968072488942, 3.0822647851901, 3.26452177629061, 
3.46223321951649, 3.66011832966054, 3.85710605543097, 4.05385887531972, 
4.83943843494744, 5.05864734149161, 5.25501778319145, 5.38941130574907, 
5.88571117751377, 6.5116611852713, 6.98632496342285, 7.21816245728101, 
7.73244825971004, 7.80401007592906, 8.34648625541999, 9.83184090479964, 
10.8324874884172, 11.3060100107816, 12.3048113953808, 13.1300123358331
), .Dim = c(99L, 2L), .Dimnames = list(NULL, c("Theta", "q(x)_(Theta)"
)))

This is my q_theta(x) function that I estimated in R. One of the question I have is:

• a> If x is a standard normal distribution this integral is zero; Right?

• b> Otherwise, in my case, the integral is not zero. How do I treat the q_1-Theta(x)? Its simply the sort(matrix[,"q(x)_(Theta)"],decreasing=TRUE) ?

And the integration would be:

sintegral(thau[1:50], (matrix[,"q(x)_(Theta)"][1:50] - sort(matrix[,"q(x)_(Theta)"],TRUE)[1:50])[1:50])$value

The median would be a comun point of this two functions. Right?

Thanks.

回答1:

Recall your previous post Building a function by defining X and Y and then Integrating in R, we build a linear interpolation function

## note `rule = 2` to enable "extrapolation";
## otherwise `rule = 1` gives `NA` outside [0.01, 0.5]
integrand <- approxfun(mat[, 1], y, rule = 2)

Then we can perform numeric integration on [0, 0.5]:

integrate(integrand, lower = 0, upper = 0.5)
# -5.594405 with absolute error < 4e-04

Now for a>, let's have a proof first.

Note, your quantile function is not for normal distribution, so this result does not hold. You can actually verify this

quant <- approxfun(mat[, 1], mat[, 2], rule = 2)
integrate(quant, lower = 0, upper = 0.5)
# -3.737973 with absolute error < 0.00029

Compared with previous integration result -5.594405, the difference is not a factor of 2.