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 не допускаются.