1

I have set some initial parameters for the moments of the gamma and lognormal distribution and applied the Kolmogorov-Smirnov test to obtain the p-value. My aim is to show a plot of p-values against N for varying N. Lets say between 5 and 1000. How would I achieve this?

mean <- 10
var <- 40
N <- 100
gamsample <- rgamma(N, shape=mean^2/var, rate=mean/var)
lnsample <- rlnorm(N, meanlog=log(mean)-log(1+mean^2/var)/2,
                          sdlog=sqrt(log(1+(mean^2/var))))

ks.test(gamsample, lnsample)$p.value

1 Answer 1

6

Use sapply():

mean <- 10
var <- 40
myP <- function(N) {
  gamsample <- rgamma(N, shape=mean^2/var, rate=mean/var)
  lnsample <- rlnorm(N, meanlog=log(mean)-log(1+mean^2/var)/2,
                     sdlog=sqrt(log(1+(mean^2/var))))

  ks.test(gamsample, lnsample)$p.value
}
N <- 5:1000
sapply(N, myP)

Or (more secure) vapply(N, myP, FUN.VALUE = 9.9).
If you want a dataframe as result you can do:

results <- data.frame(n=N, pval=sapply(N, myP))

For a plot @AndreElrico proposed in his (deleted) answer:

ggplot2::ggplot(results, aes(x=n,y=pval)) + geom_point()
5
  • Nitpicks: avoid dots in names, they mess with automatic S3 method registration. sapply is error-prone as its output type is unpredictable. Prefer vapply or lapply as appropriate. Aug 15, 2018 at 9:21
  • @KonradRudolph I changed the names.
    – jogo
    Aug 15, 2018 at 9:24
  • @jogo thank you for your answer, very educational. I have edited my question with an additional query. Not sure if I should ask another question because it's related to this but I would like to see your take on it.
    – Ali
    Aug 15, 2018 at 14:10
  • meta.stackoverflow.com/questions/275138/… ... so ask an new question and put a link to this one.
    – jogo
    Aug 15, 2018 at 14:16
  • Okay, thanks. I should have included it within this question but thought of it afterwards.
    – Ali
    Aug 15, 2018 at 14:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.