0

I am trying to use trapezoid integration, and as a result integrating the survival function(1-F(x)) of a mixture of lognormal distribution, using the code below:

for(i in 1: nrow(mydf)){
mydf$Ex[i]<- trapzfun(function(x){(1-pnorm((log(x)-mydf$mu1[i])/mydf$sd1[i]))*mydf$pmix1[i]+(1-pnorm((log(x)-mydf$mu2[i])/mydf$sd2[i]))*mydf$pmix2[i]},a=0,b=1)

mu1,sd1,pmix1 is my mean,std in logscale. pmix1 and pmix2 are mixing proportions.

While this is running well, it is taking about 3-5 hours of run time. I guess "for loops" are not a good method, and I am pretty new to R. I did try to use an apply function:

mixture<- function(x){
  trapzfun(function(x){(1-plnorm(x,mydf$mu1,mydf$sd1))*mydf$pmix1+(1-plnorm(x,mydf$mu2,mydf$sd2))*mydf$pmix2},a=0,b=1)
}

[Note the similarity with plnorm and pnorm using the transformation I applied]

mydf$Ex_1<- apply(mydf,1,mixture)

It is returning list(value = c(0.055257000747731, 0.055257000747731, 0.00.....

Can you please help me on this problem!

12b345b6b78
  • 995
  • 5
  • 16
user3807808
  • 21
  • 1
  • can you try `sapply` instead of `apply`? Looks like you are looking for an unlisted output – Vivek Kalyanarangan Dec 01 '18 at 05:47
  • Thanks Vivek! When I run the sapply on top 100 rows using: > test$Ex_1<- sapply(test,mixture) Error in `$<-.data.frame`(`*tmp*`, Ex_1, value = list(c(0.00970371416069119, : replacement has 3 rows, data has 100 – user3807808 Dec 01 '18 at 06:23
  • `sapply` is not "safe"; I suggest that you try: `unlist(lapply(seq_len(nrow(mydf)), function(i){ ...smth.... }))` or `vapply` (you also supply the return format) or `mclapply` (a parallelized version of lapply). – LostIT Dec 01 '18 at 09:12
  • Could you make your problem reproducible by sharing a sample of your data and the code you're working on so others can help (please do not use `str()`, `head()` or screenshot)? You can use the [`reprex`](https://reprex.tidyverse.org/articles/articles/magic-reprex.html) and [`datapasta`](https://cran.r-project.org/web/packages/datapasta/vignettes/how-to-datapasta.html) packages to assist you with that. See also [Help me Help you](https://speakerdeck.com/jennybc/reprex-help-me-help-you?slide=5) & [How to make a great R reproducible example?](https://stackoverflow.com/q/5963269) – Tung Dec 02 '18 at 07:53

0 Answers0