################################
### PROPAGATION INCERTITUDES ###
################################


f <- function(u) {
	v <- 2*cos(u^2+1)*atan(exp(u))
	return(v)
}

seuil <- 0.6

par(mfrow=c(1,2))
plot(t <- seq(-1,1,length=100), f(t), type="l", xlab="x", ylab="y", col="blue", lwd=2)
abline(h=-1, v=-1)
abline(h=seuil, col="red")

N <- 20

for (i in 1:N) {
	x <- runif(1, min=-1, max=1)
	y <- f(x)
	if (N<30) {
		lines(c(x,x), c(-1,y), col="black", lty="dotted")
		lines(c(-1,x), c(y,y), col="black", lty="dotted")
	}
	points(x, -1, pch=20, col="black")
	points(-1, y, pch=20, col="black")
}

N <- 10000
y <- f(runif(N, -1, 1))
hist(y, freq=FALSE, main=paste("Histogramme de y (pour ", N, " simulations)", sep="")); abline(v=seuil, col="red")
p <- sum(y>seuil)/N
cat("\nProbabilité estimée que f(X) dépasse ", seuil, " avec X uniforme sur [-1,1]: ", p, "\n\n")


#########################
### SIMULATION DISQUE ###
#########################

# simulation uniforme dans le disque D(0,R)

# methode incorrecte ...
# U : uniforme sur [0,R]
# V : uniforme sur [0,2*pi]

R <- 2
N <- 2000

S <- matrix(runif(2*N),N,2)
U <- S[,1]
V <- S[,2]*2*pi
X <- R*U*cos(V)
Y <- R*U*sin(V)

par(mfrow=c(1,3))
plot(X,Y, main="Simulation incorrecte en coordonnées polaires")
theta_aux <- seq(from=0, to=2*pi+0.1,by=0.1)
lines(R*cos(theta_aux), R*sin(theta_aux))

# methode correcte avec les coordonnées polaires
# U : loi à trouver
# V : uniforme sur [0, 2*pi]

# A VOUS DE JOUER !!

plot(0,0, main="Simulation correcte : à vous de jouer !", xlim=c(-R,R), ylim=c(-R,R))
theta_aux <- seq(from=0, to=2*pi+0.1,by=0.1)
lines(R*cos(theta_aux), R*sin(theta_aux))

# simulation par rejet
A <- runif(2*N, min=-R, max=R)
X <- A[1:N]
Y <- A[(N+1):(2*N)]
I <- (X^2 + Y^2 <= R^2)  # indices pour lesquels la condition est vérifiée
X.disque <- X[I]
Y.disque <- Y[I]
plot(X.disque, Y.disque, main="Simulation par rejet")
theta_aux <- seq(from=0, to=2*pi+0.1,by=0.1)
lines(R*cos(theta_aux), R*sin(theta_aux))