0 
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;

我正在尝试使用蒙特卡罗集成与 proc iml 来估计 x**2 - 3*x**3 + 5x 小于 1 的概率。我做错了什么？顺便说一句，Do 循环是不允许的。

4

2 回答 2

1

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
于 2014-05-09T16:14:45.920 回答
0

这可能是一个家庭作业问题，所以我不会给出完整的答案。但是，我的第一个提示是考虑 Y 应该是什么。它应该是一个指示变量（0/1），指示随机选择的观察是否落在平面的某个区域内。

我的第二个提示是记住你不能在表达式的右侧使用未定义的变量，所以 'y + (f(x)<1)' 不是有效的向量表达式。祝你好运！

于 2014-05-09T10:11:36.250 回答