CODEWITHSUNDEEP
r
Elekko
Elekko

Reputation: 109

Collect random dots under normal distribution using R

I'm trying to write an R script by first plotting standard normal distribution curve and then uniformly generating dots on the plot. I want to collect all the amount of dots that come under the standard normal PDF. But when I calculate the amount, it is not the same as displayed on the plot.

Here is my code,

N = 20
normpdf = function(x) {
    f = (1/sqrt(2*pi))*exp(x^2/(-2))
    return(f)
}        

x = seq(-5,5,length=N)
normprobdist = normpdf(x)    
plot(x,normprobdist,type="l")
par(new=TRUE) 
u = NULL
dots_undercurve = NULL
for (i in 1:length(x)) { 
     u[i] = runif(x,0,1)
     if (u[i]<=normprobdist[i]) {
                dots_undercurve[i] = TRUE
     } 
     else {
              dots_undercurve[i] = FALSE
     }

}

sum(dots_undercurve)
plot(x,u,xaxt='n', yaxt = 'n', ann=FALSE)

The sum(dots_undercurve) displays different number of dots under the curve, than seen on the plot. My programming knowledge is quite bad, so please if someone could find whats wrong with my script? Thanks

Upvotes: 0

Views: 143

Answers (2)

Philippe
Philippe

Reputation: 204

Many errors here: I thinks what you want is something like that:

rm(list=ls(all=TRUE))
options(warn=-1)
N = 20
normpdf = function(x) {
  f = (1/sqrt(2*pi))*exp(x^2/(-2))
  return(f)
}
x = seq(-5,5,length=N)
normprobdist = normpdf(x)    
par(new=TRUE) 
u = runif(N,0,1)
dots_undercurve = NULL
for (i in 1:length(x)) { 
  if (u[i]<=normprobdist[i]) {
  dots_undercurve[i] = TRUE
} 
else {
  dots_undercurve[i] = FALSE
} 
}
sum(dots_undercurve)
plot(x,u)
lines(x,normprobdist,type="l")

enter image description here

Upvotes: 0

Roland
Roland

Reputation: 132576

You don't need a for loop or to define the pdf yourself:

N <- 1e6
set.seed(42) #for reproducibility of random numbers
x <- runif(N, -5, 5)
y <- runif(N, 0, 1)

yunder <- y < dnorm(x) #which dots are under curve
sum(yunder)/N * (10 * 1) #the expected value is about 1
#[1] 0.99914

curve(dnorm, -5, 5, ylim=c(0,1))
points(x, y, col=c("red", "green")[yunder+1], pch=".")

resulting plot

Upvotes: 1

Related Questions