dgeom(X, prob = p)            
# p is the probability of success

dgeom(X, prob = p, log=TRUE)  
# log of the probability instead

dhyper(X, m, n, k)

dhyper(X, m, n, k, log = TRUE)

metropolis.hastings <- function(f,  # the target distribution

                                g,  # the proposal distribution

                                rg, # a sample from the proposal distribution

                                x0, # initial value for chain, in R it is x[1]

                                chain.size=1e5,  # chain size

                                burn.perc=0.1) { # burn in percentage

  x <- c(x0, rep(NA,chain.size-1))  # initialize chain

  for(i in 2:chain.size)   {

    y <- rg(x[i-1])                 # generate Y from g(.|xt) using sampler rg

    alpha <- min(1, f(y)*g(x[i-1],y)/(f(x[i-1])*g(y,x[i-1])))

    x[i] <- x[i-1] + (y-x[i-1])*(runif(1)<alpha)  # update step

  }

  # remove initial part of the chain before output result

  x[(burn.perc*chain.size) : chain.size] 

}

a<-2.7; b<-6.3; size<-1e4

f  <- function(x)   dbeta(x,a,b)
rg <- function(x)   runif(1,0,1)
g  <- function(x,y) 1 # i.e., dunif(x,0,1)

X <- metropolis.hastings(f,g,rg,x0=runif(1,0,1),chain.size=size)

par(mfrow=c(1,2),mar=c(2,2,1,1))
hist(X,breaks=50,col="blue",main="Metropolis-Hastings",freq=FALSE)
curve(dbeta(x,a,b),col="sienna",lwd=2,add=TRUE)
hist(rbeta(size,a,b),breaks=50,col="grey",main="Direct Sampling",freq=FALSE)
curve(dbeta(x,a,b),col="sienna",lwd=2,add=TRUE)

library(VGAM)

sigma <- 4
f  <- function(x)   drayleigh(x,sigma) # the target distribution
rg <- function(x)   rchisq(1,df=x)     # a sample from proposal g(.|x)
g  <- function(x,y) dchisq(x,df=y)     # the pdf at g(x|y)

X <- metropolis.hastings(f,g,rg,x0=rchisq(1,df=1),chain.size=size)

par(mfrow=c(1,2),mar=c(2,2,1,1))
hist(X,breaks=50,col="blue",main="Metropolis-Hastings",freq=FALSE)
curve(drayleigh(x,sigma),col="sienna",lwd=2,add=TRUE)
hist(rrayleigh(size,sigma),breaks=50,col="grey",main="Direct Sampling",freq=FALSE)
curve(drayleigh(x,sigma),col="sienna",lwd=2,add=TRUE)

          "@R_Experts"

meshgrid <- function(a, b) {

  list(x = outer(b * 0, a, FUN = "+"), 

       y = outer(b, a * 0,FUN = "+"))

}

#install.packages("rgl")

library(rgl)

c = 0.2

b = 1.7

theta = seq(0, 2 * pi, length = 40 * 4)

phi = seq(0, pi, length = 40 * 4)

theta1 = meshgrid(theta, phi)$x

phi2 = meshgrid(theta, phi)$y

x = (1 + c * phi2) * sin(phi2) * cos(theta1)

y = (1 + c * phi2) * sin(phi2) * sin(theta1)

z <- b * cos(phi2)

surface3d(x, y, z, color = rainbow(10))

par3d(zoom = 0.7)

p

q

d

r