How to run a monte carlo simulation from a custom distribution in R

Question

I would like to pull 1000 samples from a custom distribution in R

I have the following custom distribution

library(gamlss)
mu    <- 1    
sigma <- 2 
tau   <- 3   
kappa <- 3
rate  <- 1
Rmax  <- 20

x <- seq(1, 2e1, 0.01)
points <- Rmax * dexGAUS(x, mu = mu, sigma = sigma, nu = tau) * pgamma(x, shape = kappa, rate = rate)
plot(points ~ x)

How can I randomly try Monte Carlo simulations?

My first attempt was the following code, which created a histogram shape that I did not expect.

hist(sample(points, 1000), breaks = 51)

This is not what I was looking for, as it does not match the same distribution as in pdf format.

Answer 1

If you need a Monte Carlo simulation, you will need to select a large number of times from the distribution, rather than taking a large sample once.

points , 400, , . plot(points ~ x). , points . , . , x y plot(points ~ x). points , , 1000 1900 . points ( ):

hist(points, 100)

100 , .

, , , ( x). , points 2 , 1. , , plot(points ~ x) , 2, 0.5 1.5. plot(points ~ x). , ( , ) 0 0.25. , , - , :)

- , - :

samples <- replicate(1000, sample(points, 100, replace = TRUE))

, points ,

Answer 2

ECDF :

 ecd.points <- ecdf(points)
 invecdfpts <- with( environment(ecd.points), approxfun(y,x) )
 samp.inv.ecd <- function(n=100) invecdfpts( runif(n) )
 plot(density (samp.inv.ecd(100) ) )
 plot(density(points) )
 png(); layout(matrix(1:2,1)); plot(density (samp.inv.ecd(100) ),main="The Sample" )
  plot(density(points) , main="The Original"); dev.off()

Answer 3

( ) :

library(gamlss)
fun <- function(x, mu = 1, sigma = 2, tau = 3, kappa = 3, rate = 1, Rmax = 20)
  Rmax * dexGAUS(x, mu = mu, sigma = sigma, nu = tau) * 
  pgamma(x, shape = kappa, rate = rate)

MCMC ( Monte Carlo). ,

simMCMC <- function(N, init, fun, ...) {
  out <- numeric(N)
  out[1] <- init
  for(i in 2:N) {
    pr <- out[i - 1] + rnorm(1, ...)
    r <- fun(pr) / fun(out[i - 1])
    out[i] <- ifelse(runif(1) < r, pr, out[i - 1])
  }
  out
}

It starts with a dot initand gives a Ndraw. This approach can be improved in many ways, but I'm just going to start the form init = 5, turn on the burn period of 20,000 and select every second draw to reduce the number of repetitions:

d <- tail(simMCMC(20000 + 2000, init = 5, fun = fun), 2000)[c(TRUE, FALSE)]
plot(density(d))

Answer 4

Here is another way to do this, which is extracted from R: Generate data from a probability density distribution and How to create a distribution function in R? :

x <- seq(1, 2e1, 0.01)
points <- 20*dexGAUS(x,mu=1,sigma=2,nu=3)*pgamma(x,shape=3,rate=1)
f <- function (x) (20*dexGAUS(x,mu=1,sigma=2,nu=3)*pgamma(x,shape=3,rate=1))
C <- integrate(f,-Inf,Inf)

> C$value
[1] 11.50361

# normalize by C$value
f <- function (x) 
(20*dexGAUS(x,mu=1,sigma=2,nu=3)*pgamma(x,shape=3,rate=1)/11.50361)

random.points <- approx(cumsum(pdf$y)/sum(pdf$y),pdf$x,runif(10000))$y
hist(random.points,1000)

hist((random.points*40),1000) will get scaling like your original feature.