#------------------------------------------------------------- # alcune funzioni di proposta # (anche esercizio su trasformazioni di variabili) # # GdA febbraio 2011 #------------------------------------------------------------- #-------------------------------------------------------------- # uniforme in phi e in distanza ------------------------------- proponi.unif <- function(xy.old, d) { phi <- runif(1)*2*pi r <- runif(1)*d dx <- r * cos(phi) dy <- r * sin(phi) return ( c(xy.old[1]+dx, xy.old[2]+dy) ) } # test n = 10000 x=y=rep(0,n) for (i in 1:n) { xy <- proponi.unif( c(0,0), 1) x[i] <- xy[1] y[i] <- xy[2] } # vediamo il risultato plot(x,y,cex=0.15) hist(x,nc=100) # => uniforme in phi e in d non significa uniforme in x e y!! # cosa succede in phi e d? phi = atan2(y, x) d = sqrt(x^2+y^2) plot(phi,d,cex=0.15) hist(phi,nc=100) hist(d,nc=100) #--------------------------------------------------------------- # uniforme su un quadrato -------------------------------------- proponi.unif.xy <- function(xy.old, d) { dx <- (runif(1)-0.5)*d dy <- (runif(1)-0.5)*d return ( c(xy.old[1]+dx, xy.old[2]+dy) ) } # test n = 10000 x=y=rep(0,n) for (i in 1:n) { xy <- proponi.unif.xy( c(0,0), 1) x[i] <- xy[1] y[i] <- xy[2] } # vediamo il risultato plot(x,y,cex=0.15) hist(x,nc=100) # cosa succede in phi e d? phi = atan2(x,y) d = sqrt(x^2+y^2) plot(phi,d,cex=0.15) hist(phi,nc=100) hist(d,nc=100) #--------------------------------------------------------------- # uniforme nel cerchio ----------------------------------------- proponi.unif.cerchio <- function(xy.old, d) { phi <- runif(1)*2*pi r <- sqrt(2*runif(1)) # pensare al motivo di questa trasformazione dx <- r * cos(phi) dy <- r * sin(phi) return ( c(xy.old[1]+dx, xy.old[2]+dy) ) } # test n = 10000 x=y=rep(0,n) for (i in 1:n) { xy <- proponi.unif.cerchio( c(0,0), 1) x[i] <- xy[1] y[i] <- xy[2] } # vediamo il risultato plot(x,y,cex=0.15) hist(x,nc=100) # cosa succede in phi e d? phi = atan2(y, x) d = sqrt(x^2+y^2) plot(phi,d,cex=0.15) hist(phi,nc=100) hist(d,nc=100) #------------------------------------------------------------- # gaussiana in x e in y proponi.norm <- function(xy.old, sigma) { dx <- rnorm(1, 0, sigma) dy <- rnorm(1, 0, sigma) return ( c(xy.old[1]+dx, xy.old[2]+dy) ) } # test n = 10000 x=y=rep(0,n) for (i in 1:n) { xy <- proponi.norm( c(0,0), 1) x[i] <- xy[1] y[i] <- xy[2] } # vediamo il risultato plot(x,y,cex=0.15) hist(x,nc=100) # cosa succede in phi e d? phi = atan2(y, x) d = sqrt(x^2+y^2) plot(phi,d,cex=0.15) hist(phi,nc=100) hist(d,nc=100) # cosa ricorda?