#-------------------------------------------------------------
# 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?

Created by Pretty R at inside-R.org