[R] integration in R?

[Martyn Plummer]

On 20-Sep-00 Bill Simpson wrote:
> I have the following problem.
> 
>    inf
> p= Int  [dnorm(x-d) * pnorm(x) ^ (m-1)]  dx
>   -inf
> 
> Given p and m, I want to find d. For example,
> 
> p     m       d
> .9    2       1.81
> 
> Is there a way to solve this problem in R? Ideally I would have a
> function with arguments p and m, and output d.

This function returns an approximation to the integral, which
should work well.

foo <- function(d, mm, N=10000, llim=-5, ulim=5)
{
   x <- seq(from=llim, to=ulim, length=N)
   mean(dnorm(x) * pnorm(x+d)^(mm-1)) * (ulim - llim)
}

To solve the inverse problem, create another function that returns
zero at the solution to your equation

bar <- function(d, p, mm, N=10000, llim=-5, ulim=5) {
   x <- seq(from=llim, to=ulim, length=N)
   p - mean(dnorm(x) * pnorm(x+d)^(mm-1)) * (ulim - llim)
}

and then feed it into uniroot()

R> uniroot(bar, interval = c(0,5), p=0.9, mm=2)
$root
[1] 1.81312

$f.root
[1] -5.655575e-07

$iter
[1] 7

$estim.prec
[1] 6.103516e-05

I had to rename the argument "m" to "mm" to get uniroot() to accept this
function. Without this, the "m" argument is skimmed off as the "maxiter"
parameter of uniroot(), and is not passed on to bar().

Martyn
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._