Использование proc iml для интеграции по методу Монте-Карло

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


person user3547583    schedule 09.05.2014    source источник
comment
Я предполагаю, что более важный вопрос, где будет, как мне сделать (x2 - 3*x3 + 5x ‹ 1) счет в матрице?   -  person user3547583    schedule 09.05.2014


Ответы (2)


Здесь много чего не так. Как говорит Рик, это похоже на домашнюю задачу, поэтому я не буду публиковать решение. Вот некоторые вещи, о которых стоит подумать.

  1. Обычное судебное разбирательство, вероятно, не то, что вы ищете. Ваша функция определена на всей числовой прямой. Получение x~N(0,1), где abs(x) > 3, маловероятно.
  2. Рик упоминает, что вы используете y перед его определением.
  3. Посмотрите на поэлементные операции вместо матричных операций. Используйте их в своей функции.
  4. 5x недопустимая операция.
  5. Рассмотрите возможность сокращения множества [:], чтобы получить среднее значение.
  6. Я не думаю, что ваше уравнение для SE правильно. http://en.wikipedia.org/wiki/Standard_error
person DomPazz    schedule 09.05.2014

Вероятно, это домашнее задание, поэтому я не буду давать полный ответ. Тем не менее, мой первый совет — подумать о том, каким должно быть Y. Предполагается, что это индикаторная переменная (0/1), указывающая, попадает ли случайно выбранное наблюдение в определенную область плоскости.

Мой второй совет: помните, что вы не можете использовать неопределенную переменную в правой части выражения, поэтому 'y + (f(x)‹1)' не является допустимым векторным выражением. Удачи!

person Rick    schedule 09.05.2014