Does this do any better, or perhaps offer other ideas? I could not understand what the "Hct" and step()'s were doing from the description of your desired result, so they were omitted from geq(). I also could not understand why you were rounding off the integrals or why there wee two integral terms. I think the problem may be (in part) that integrate() returns a list unless you specify the value column. When I tested your function or my earliest efforts at fixing it I continued to get a single value despite passing it a sequence or a matrix.
surviv<-function(x){k=.001;alpha=0.3;tau=70; exp(-k*x)/(1+exp(alpha*(x-tau)))} geq<-function(n){ integrate(surviv, 0,n)$value } Nmax=250 gxt<-matrix(nrow=Nmax) gxt<-as.matrix(1:250) # there may be better methods o making this matrix geq.results<-geq.results<-apply(gxt,1,geq) geq.results plot(gxt, geq.results) I suspect that I have not understood what you wanted from geq(.) since the plot looks rather "boring". I am guessing that Hct0 is a hematocrit value. And geq basically returns a value very close to its input up until tau and then stays at approximately tau for the rest of time. -- David Winsemius, MD Robert Kalicki wrote: > Hello everybody! > > I have problems with integrating my function. > > My primary function is a “survival function” of the following type: > > surviv <- exp(-k*x)/(1+exp(alpha*(x-tau))) > > I would like to integrate this function over a defined range and obtain a > vector with all the values from integrate(surviv, 0, x) for x <- 1:N. > > I unfortunately obtain just a scalar value corresponding to the last > integrated x. > > How do I have to proceed if I want to obtain a vector with all the results? > > I have tried to solve this problem by adding a “while loop”. It works but > the expression becomes to complex if I want to perform further optimization. > > Have any one a solution for this problem? > > Thanks a lot in advance. > > > > Example > > # General equation > > beta1 <- 0.001 > beta2 <- 0 > Hct0 <- 0.27 > tau <- 70 > k <- 0.001 > alpha <- 0.3 > T2 <-40 > > step <- function(x){ > if(x>=0) step <- 1 > else step <- 0 > > } > Nmax <- 250 > n <- 1 > while(n<Nmax){ > surviv <- function(n){ > surviv <- exp(-k*n)/(1+exp(alpha*(n-tau))) > } > geq <- function(n){ > geq <- > Hct0+beta1*as.integer(integrate(surviv,0,n)[1])+step(n-T2)*(beta2-beta1)*as. > integer(integrate(surviv,0,n-T2)[1]) > } > plot(n, geq(n), type="p", pch=19, xlim=c(0,250), > ylim=c(0.25,0.35), xlab="time in days", ylab="Hct") > > par(new=T) > print(c(n, surviv(n), geq(n), > as.integer(integrate(surviv,0,n)[1])) > n <- n + 1 > } > > Robert M. Kalicki, MD > > Department of Nephrology and Hypertension ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.