R : приближение вероятности

Пусть X и Y — две независимые случайные величины с функциями плотности f_x(x) = x*exp(0,5*x^2) где x>0 и f_y(y) = 0,5 с y в [-1,1] .

Я пытаюсь аппроксимировать вероятность P(Z>1) где Z =X+Y ; вот мой код R:

u <-runif(10000, min=0, max=1000000); 
v<-runif(10000, min=-1, max=1);
f<-function(x){x*exp(-0.5*x^2)} ; 
g<-function(x){1/2} ; 
for ( i in 1:10000) {
I[i] = integrate ( f ,lower = 0 , upper = u[i]) ; 
J[i] = integrate ( g, lower =-1 , upper = v[i]) ; 
}
mean( (I+J)>1 ) ;

Я получаю эту ошибку:

Error in integrate(g, lower = -1, upper = v[i]) : 
  evaluation of function gave a result of wrong length
In addition: Warning message:
In I[i] = integrate(f, lower = 0, upper = u[i]) :
  number of items to replace is not a multiple of replacement length

r approximation probability
person aouakki    schedule 23.12.2015    source источник
comment
Это может оказаться полезным.   -  person ytk    schedule 23.12.2015
comment
отрицательный знак в exp опечатка? Поскольку вы упоминаете в описании, что ваша функция плотности: f_x (x) = x * exp (0,5 * x ^ 2)   -  person astrosyam    schedule 23.12.2015
comment
Может быть, описание неправильное, я думаю   -  person astrosyam    schedule 23.12.2015
comment
Вы понимаете, что возвращаемое значение из integrate - это не число, а список с именованными элементами????   -  person IRTFM    schedule 23.12.2015
comment
@ mra68 В приведенном выше коде есть (как минимум) две проблемы. A) результат от integrate не будет числовым, поэтому вы не можете его ни к чему добавить, и B) вы оба не читаете страницу справки, где четко указаны требования функции, переданной integrate. Мой (теперь перевернутый) голос против вашего ответа заключается в том, что я считаю бесполезным обвинять программное обеспечение в наличии ошибок просто потому, что вы не знаете, как его использовать.   -  person IRTFM    schedule 23.12.2015

Ответы (1)


arrow_upward
0
arrow_downward

Замените g <- function(x){1/2} на g <- Vectorize(function(x){1/2} ) и integrate(...) на integrate(...)$value:

u <-runif(10000, min=0, max=1000000)
v<-runif(10000, min=-1, max=1)
f<-function(x){x*exp(-0.5*x^2)} 
g<-Vectorize(function(x){1/2 + 1e-12*x})

I <- rep(NA,10000)
J <- rep(NA,10000)

for ( i in 1:10000) {
  I[i] = integrate ( f ,lower = 0 , upper = u[i])$value
  J[i] = integrate ( g , lower =-1 , upper = v[i])$value 
}
mean( (I+J)>1 )

.

> mean( (I+J)>1 ) 
[1] 0.0048
>
person mra68    schedule 23.12.2015
comment
Это не ошибка в integrate. Вы дали ему функцию, которая возвращает только один элемент. Прочтите страницу справки. - person IRTFM; 23.12.2015
comment
Это также было бы успешным и, вероятно, более эффективным: g<-function(x){rep(1/2, times=length(x))} - person IRTFM; 23.12.2015