CODEWITHSUNDEEP
sasmontecarlosas-iml
user3547583
user3547583

Reputation: 1

Using proc iml to do monte carlo integration

proc iml;

call randseed(4545); * initialize the stream (like streaminit);
x = J(5000,1,.); * pre-allocate space for random numbers;
call randgen(x,'normal',0,1); * fill x with N(0,1) deviates;
y = y + (x**2 - 3*x**3 + 5x < 1);
p = y / 5000;  * MEAN acts on each matrix column;
se = sqrt(y*(1-y)/5000); * VAR, but not STD or STDERR, also acts on columns;
print "IML STEP: estimated probability is" p
"with standard error" se;  * use PRINT, not PUT;

I'm trying to use monte carlo integration with proc iml to estimate the probability that x**2 - 3*x**3 + 5x is less than 1. What am I doing wrong? Do loops are not allowed by the way.

Upvotes: 0

Views: 386

Answers (2)

DomPazz
DomPazz

Reputation: 12465

Lots of things wrong here. AS Rick says, this has the feel of a homework problem, so I won't post the solution. Here are some things to think about.

  1. A normal trial is probably not what you are looking for. Your function is defined on the entire number line. Getting a x~N(0,1) where abs(x) > 3 is not likely.
  2. Rick mentions that you are using y before defining it.
  3. Look at element wise operations instead of matrix operations. Use these in your function.
  4. 5x is not a valid operation.
  5. Consider set reduction [:] to do the mean.
  6. I don't think your equation for the SE is correct. http://en.wikipedia.org/wiki/Standard_error

Upvotes: 1

Rick
Rick

Reputation: 1210

This is probably a homework problem, so I am not going to give the complete answer. However, my first hint is to think about what Y is supposed to be. It is supposed to be an indicator variable (0/1) that indicates whether a randomly chosen observation falls within a certain region of the plane.

My second hint is to remember that you cannot use an undefined variable on the right hand side of an expression, so 'y + (f(x)<1)' is not a valid vector expression. Good luck!

Upvotes: 0

Related Questions