[R] How to find when a value is reached given a function?

Abby Spurdle @purd|e@@ @end|ng |rom gm@||@com
Tue Jan 26 01:52:43 CET 2021


You could use a spline to interpolate the points.
(And I'd consider increasing the number of points if possible, say to 200).

Then use a root finder, such as uniroot(), to solve for
f(i) - k
Where, k (a constant), would be 1e6, based on your example.

There are a number of variations on this approach.
My kubik package provides a solve method, and can impose some constraints.

----
library (kubik)
f <- chs (1:45, round (PCR_array),
    constraints = chs.constraints (increasing=TRUE) )
plot (f)

sol <- solve (f, 1e6)
abline (v=sol, lty=2)
sol
----

Note that I had to round the values, in order to impose a
non-decreasing constraint.

Also note that I've just used the 45 points.
But re-iterating, you should increase the number of points, if possible.


On Mon, Jan 25, 2021 at 8:58 AM Luigi Marongiu <marongiu.luigi using gmail.com> wrote:
>
> Hello
> I am trying to simulate a PCR by running a logistic equation. So I set
> the function:
> ```
> PCR <- function(initCopy, dupRate, Carry) {
>   ROI_T = initCopy
>   A = array()
>   for (i in 1:45) {
>     ROI_TplusOne <- ROI_T * dupRate * (1 - ROI_T/Carry)
>     A[i] <- ROI_TplusOne
>     ROI_T <- ROI_TplusOne
>   }
>   return(A)
> }
> ```
> Which returns an array that follows the logistic shape, for instance
> ```
> d <- 2
> K <- 10^13
> A_0 <- 10000
> PCR_array <- PCR(A_0, d, K)
> plot(PCR_array)
> ```
> Given the formula `ROI_TplusOne <- ROI_T * dupRate * (1 -
> ROI_T/Carry)`, is it possible to determine at what time point `i` a
> given threshold is reached? For instance, what fractional value of i
> returns 1000 000 copies?
> Thank you
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.



More information about the R-help mailing list