Tengo dos variables aleatorias $s\sim \mathcal{N}(\nu,\sigma)$ y $a\sim \mathcal{U}(0,A)$, $0<A<1$, calcular una tercera r.v. $t=(1-a)/s$ y quiere encontrar su distribución $p(t|\nu,\sigma,A)$.
Mi razonamiento es, que para cualquier valor de $a$, no es exactamente una $s=(1-a)/t$ de manera tal que el $(a,s)$ par producirá $t$, y por lo tanto tendrá que integrar sobre las probabilidades para estos valores:
$$p(t|\nu,\sigma,A) = \int_0^A p_a(x)p_s((1-x)/t)\,\mathrm dx\\$$
que puede ser expresado como la suma de dos funciones de error.
Pero ya cuando el control de este paso numéricamente, tengo una discrepancia entre la estima (por muestreo) y calculado de distribución:
¿Por qué el histograma y el calculado de distribución no coinciden?
Este es el código que he usado:
n=10000000
A=0.7
nu=.7
sigma=.1
# sampling from the target distribution
s=rnorm( n, mean=nu, sd=sigma )
a=runif( n, min=0, max=A )
t=(1-a)/s
hist(t,200, freq=F, xlim=c(0,5), ylim=c(0,2.5))
# plot the analytical result
analytic <- function(t, A, nu, sigma ){
tmp <- function(x, A, t, nu, sigma){
return (1/A*dnorm( (1-x)/t, mean=nu, sd=sigma) )
}
return( integrate( tmp, 0, A, A=A, t=t, nu=nu, sigma=sigma)$value)
}
x=seq(0,5,by=.01)
y=rep(0,length(x))
for( i in seq(1,length(x)) ){
y[i]=analytic(x[i], A, nu, sigma)
}
lines( x, y, col="red", type="l")
