kjør simulering for å finne sannsynligheter (OBIG II, oppgave 2)

Svært kjekt med eksakte svar, på de to sannsynligheter jeg spør om, i oppgave 2, via flinke integraler. Men man kan også ty til simulering -- kjør dette (i R), jeg lager en haug med (x, y) fra oppgitt fordeling, og kan "lese ut" enhver sannsynlighet. 
Nils Lid Hjort

# (X,Y) uniform on circle,
sim <- 10^4
xsim0 <- runif(sim,-1,1)
ysim0 <- runif(sim,-1,1)

keep <- 0*(1:sim)
for (ss in 1:sim)
{ keep[ss] <- 1*(xsim0[ss]^2 + ysim0[ss]^2 <= 1) }

mean(keep) # should be pi / 4
print( round(c(mean(keep),pi/4),4) ) # good

xsim <- xsim0[keep==1]
ysim <- ysim0[keep==1]
plot(xsim,ysim)

phat <- mean(1*(xsim <= 0.50)) # close to 0.804
qhat <- mean(1*(xsim >= 0.50)*(ysim >= 0.50))

print(round(c(phat,qhat),4)) # close to the Real Things if sim is big

# with sim = 10^6: 0.8044 and 0.0248; exact: 0.8045 and 0.0251