0

I am using R. I have made a random function (using monte carlo integration) that approximates the arctan function (equal to the integral of 1/(1+x^2)). It seems to work well and gives accurate approximations. For example 3*arctan(sqrt(3)) should give pi, and below I get a close approximation.

> f=function(x)
+ {y=runif(10^6,0,x); return(x*mean((1/(1+y^2))))}
> f(sqrt(3))*3
[1] 3.140515

However, when I use the function on a sequence of numbers, the answers seem to come out wrong. Below the correct values are given by the atan function, but the f function gives the wrong values when used on the same vector.

> atan(c(1,1.5,2))
[1] 0.7853982 0.9827937 1.1071487
> f(c(1,1.5,2))
[1] 0.6648275 0.9972412 1.3296550
> f(1)
[1] 0.7852855
> f(1.5)
[1] 0.9828134
> f(2)
[1] 1.107888

Notice how when I use f on 1, 1.5 and 2 individually, it works, but not in vector form? What's going on here? I have tried running f a few times on the vector and the wrong values are pretty consistent to 2 decimal places.

user124249
  • 79
  • 5
  • 2
    Using `mean()` inside your function interferes with vectorization. You can wrap your function in `Vectorize()` to fix it. `f=Vectorize(function(x) { ...})` – MrFlick Jan 15 '16 at 23:12
  • @MrFlick Thanks, that fixed it. If you comment this as an answer, I could mark this question as solved. – user124249 Jan 15 '16 at 23:18
  • Vectorize() is basically a wrapper for mapply(), IIRC. If your function takes one vector and should return vector, you could get it cheaper with sapply() – Severin Pappadeux Jan 16 '16 at 00:43

1 Answers1

3

Using mean() inside your function interferes with vectorization. You can wrap your function in Vectorize() to fix it. For example

f <- Vectorize(function(x) {y=runif(10^6,0,x); return(x*mean((1/(1+y^2))))})
MrFlick
  • 195,160
  • 17
  • 277
  • 295