-1

Trying to find integral (o to pi/2) of x^2 cosx using montecarlo method. This is my first time so need some direction. To generate random number should I transform the limit of the integral to (0,1) or can I generate random numbers with the given limit (0 to pi/2) ?

something like this ?

transform the integral from (o to pi/2) to (0to 1) which tranforms the function to 1/x^2 sinx generate random number rnorm(10000,0,1)

or Is there a way to generate random number like this rnorm(10000,0,1)*pi/2 with out having to transform the limit of the integral

greenH0rn
  • 65
  • 2
  • 9
  • People love to downvote here (i upvote)...but show some code first. – Fernando Oct 13 '13 at 22:05
  • Can the down voters please comment why so I can use the feed back in future ? – greenH0rn Oct 13 '13 at 22:42
  • @greenH0rn About downvotes. Your question is basically asking: can I generate uniformly distributed random numbers in given interval in R? The answer is - yes you can and it is easily findable. – sashkello Oct 14 '13 at 06:29

2 Answers2

5

You can generate random numbers uniformly in any interval you want, for example, runif(1000,0,pi/2) will generate a sample of size 1000 uniformly distributed on [0,π/2]. You definitely don't want to use rnorm here -- rnorm generates normally distributed data, not uniformly distributed data.

You could do your monte carlo simulation like this:

> f<-function(x) x^2 * cos(x)
> mean(f(runif(100000,0,pi/2)))*(pi/2)
[1] 0.4672985

Or, you can let R do the integration using integrate:

> integrate(f,0,pi/2)
0.4674011 with absolute error < 5.2e-15
mrip
  • 14,913
  • 4
  • 40
  • 58
  • 1
    Because you need to sample the entire space uniformly, there's no point in calculating the same point twice in this case (see http://en.wikipedia.org/wiki/Mean_value_theorem) – Fernando Oct 13 '13 at 22:29
  • 2
    Normally distributed data will give some parts of the interval more weight than others. Uniformly distributed data will give equal weight, which is what you want when computing an integral. – mrip Oct 13 '13 at 22:33
  • maybe a dumb question but why is not advisable to use normally distributed data ? Why is uniform distribution preferred ? also why do we need to multiply mean with pi/2 again ? – greenH0rn Oct 13 '13 at 22:40
  • When you do an ordinary integral, by definition that assumes equal weighting (i.e. uniform distribution). The normal distribution will give you a different weighting, and therefore a different (wrong) answer. You multiply the mean by pi/2 because that is the length of the interval you are integrating over. Also, these kinds of questions about integration are more appropriate for a math forum. This board is about programming questions. – mrip Oct 13 '13 at 22:43
0

My first shoot:

mc.integral = function(FUN, n.iter = 1000, interval){

  # take a sample using 'runif'
  x = runif(n.iter, interval[1], interval[2])

  # apply the user-defined function
  y = FUN(x)

  # calculate
  mean(y)*(interval[2] - interval[1])
}

example

FUN = function(x){x^2 * cos(x)}
integ = mc.integral(FUN, interval = c(0, pi/2))
print(integ)

[1] 0.4693398
Fernando
  • 7,785
  • 6
  • 49
  • 81