R ggplot2 adding function to a histogram

Question

I am trying to add this function curve to a histogram. Individually, they work. But when I try to put them on the same graph, the function messes up... I can't seem to figure out how to get them together.

# make dataframe for ggplot (can use random numbers from 0 to 10 to simulate x)
c= data.frame(x= x, cx= c(seq(from= 0.001, to= 10, by= 0.001)))

x and cx have the same number of data points.

# function for curve (alpha and beta are just constants set at 0.5)
fx= function(x){
    (beta^alpha)/((x+beta)^(alpha+1)*gamma(alpha))
}

When either the geom_histogram or the stat_function is commented out, the graph works correctly.

# graph code
h_x= ggplot(data= NULL) +
    geom_histogram(data= c, aes(x= x, y= ..density..), binwidth= 0.2, col= "purple", fill= "yellow") +
    stat_function(fun= fx, data= c, aes(x= cx)) +
    coord_cartesian(xlim= c(0:10)) +
    labs(title= "Figure 03", x= "x")
plot(h_x)

Curve by itself

curve by itself ;

Histogram and curve together

histogram and curve together

To make your example work I used set.seed(47); cc = data.frame(x= runif(10000, 0, 10), cx = c(seq(from= 0.001, to= 10, by= 0.001))) and fx = function(x, alpha = 0.5, beta = 0.5, gamma = 0.5){ (beta^alpha)/((x+beta)^(alpha+1)*gamma(alpha)) }. (I renamed the data cc because c is already the name of the most common function). With those changes to make it run, the graph looks fine. If you still have problems, please edit the code so that the problems can be reproduced. For now I am voting to close as "typo/not reproducible". — Gregor Thomas
– Gregor Thomas, Commented Aug 17, 2017 at 20:39
My guess is that somewhere you changed one (at least) of your definitions of alpha, beta, or gamma. This is why it is good practice to treat function parameters as arguments rather than having them go out looking for values in the global environment. — Gregor Thomas
– Gregor Thomas, Commented Aug 17, 2017 at 20:42

Marco Sandri · Accepted Answer · 2017-08-17 20:47:59Z

1

Like @Gregor, I made some changes to your code and the graph looks OK.
I hope it can help you.

set.seed(1)
x <- rgamma(10000,1)
df1 <- data.frame(x= x, cx= c(seq(from= 0.001, to= 10, by= 0.001)))

beta <- alpha <- 0.5
fx <- function(x) {
    print(str(x))
    (beta^alpha)/((x+beta)^(alpha+1)*gamma(alpha))
}

# graph code
h_x <- ggplot(data=df1) +
    geom_histogram(aes(x= x, y= ..density..), binwidth= 0.2, col= "purple", fill= "yellow") +
    stat_function(fun=fx, aes(x=cx), lwd=1) +
    coord_cartesian(xlim= c(0:10)) +
    labs(title= "Figure 03", x="x")
plot(h_x)

answered Aug 17, 2017 at 20:47

Marco Sandri

24.5k7 gold badges60 silver badges64 bronze badges

Sign up to request clarification or add additional context in comments.

Comments

user8480288 · Accepted Answer · 2017-08-19 15:20:36Z

0

Thanks for the help! I ended up figuring out the problem... It's because there was some large values (greater than 100) for my x values, when I removed these points, the graph looked much better!

But now my graph looks like this:

n= 10000
 i= 1
 alpha= 0.5
 beta= 0.5
 x= matrix(data= 5, nrow= n)
 lambda= matrix(data= 1.5, nrow= n)

 while (i < n) {
   x[i+1]= rexp(1, rate= lambda[i])
   lambda[i+1]= (x[i+1]+beta)^(alpha+1)*(lambda[i]^alpha)*exp(-lambda[i]*(x[i+1]+beta))

if ((lambda[i+1] < 0.00001) || (lambda[i+1] > 10)) {
  while ((lambda[i+1] < 0.00001) || (lambda[i+1] > 10)) {
    x[i+1]= rexp(1, rate= lambda[i])
    lambda[i+1]= (x[i+1]+beta)^(alpha+1)*(lambda[i]^alpha)*exp(-lambda[i]*(x[i+1]+beta))
  }
}
i= i+1
}

# data frame:
df4= data.frame(x= x[x<100], cx= c(seq(from= 0.011, to= 10, by= 0.001)))  

# graph (same function (fx) from first post):
h_x= ggplot(data= df4) +
    geom_histogram(aes(x= x, y= ..density..), binwidth= 0.2, col= "purple", fill= "yellow") +
    stat_function(fun= fx) +
    coord_cartesian(xlim= c(0:10)) +
    labs(title= "Figure 03", x= "x")
  plot(h_x)

Is there any way to make it a smooth curve? I tried scale_x_continuous but to no avail...

answered Aug 19, 2017 at 15:20

user8480288

211 bronze badge

1 Comment

user8480288 Over a year ago

Nevermind, I figured it out again.... I needed to add the argument 'n=10000' to my stat_function!

Collectives™ on Stack Overflow

R ggplot2 adding function to a histogram

2 Answers 2

Comments

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related