Rory Ellis
Rory Ellis

Reputation: 25

Simulating data for a Gompertz curve

I have a set of data that I have collected which consists of a time series, where each y-value is found by taking the mean of 30 samples of grape cluster weight.

The growth follows a Gompertz curve with formula y = a*exp(-exp(-(x-x0)/b)), with

The data:

x = c(0, 28, 36, 42, 50, 58, 63, 71, 79, 85, 92, 99, 106, 112)
y = c(0, 15, 35, 55, 62, 74, 80, 96, 127, 120, 146, 160, 177, 165).

I want to simulate more data from this, with the same number of x and y values, so that I can carry out some Bayesian analysis to find the posterior distribution of the data.

Effectively what I need is:

If there is some skeleton code where it is possible to change around the parameters, then this could potentially be very helpful for me too.

Thanks

Upvotes: 0

Views: 1395

Answers (1)

Lstat
Lstat

Reputation: 1460

Let's inspect your data

x <- c(0, 28, 36, 42, 50, 58, 63, 71, 79, 85, 92, 99, 106, 112)
y <- c(0, 15, 35, 55, 62, 74, 80, 96, 127, 120, 146, 160, 177, 165)

and fitted Gompertz curve

gFun <- function(x){
 a <- 88.8
 b <- 11.7
 x0 <- 15.1
 est <- a*exp(-exp(-(x-x0)/b))
  return(est)
}

by visualisation:

library(ggplot2)
 ggplot(ggData, aes(x=x, y=y) ) +
  geom_point() + 
  stat_function(fun=gFun, colour="blue") + 
  theme_bw()

enter image description here

This doesn't look as a good fit. However, simulating data y|x at fixed x as in the vector above can be done by adding error term. I've used normal distribution with sd=4 for illustration.

nSim <- 10

simData <- data.frame(x=c(0, rep(x[-1], each=nSim)) ) # x[-1] removes 0 from simulation
simData$y <- gFun(simData$x) + rnorm(n=nrow(simData), sd=4)

ggplot(simData, aes(x=x, y=y) ) +
 geom_point(alpha=0.4) + 
 stat_function(fun=gFun, colour="blue") + 
 scale_x_continuous(limits=c(0, max(x)) ) +
 theme_bw()

enter image description here

Upvotes: 1

Related Questions