R help - Jan 2021 - How to find when a value is reached given a function?

On 24/01/2021 2:57 p.m., Luigi Marongiu 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?
There are two answers:

The brute force answer is just to try it and count how far you need to 
go.  This is really simple, but really inefficient.

The faster and more elegant way is to solve the recursive relation for 
an explicit solution.  You've got a quadratic recurrence relation; 
there's no general solution to those, but there are solutions in special 
cases.  See https://math.stackexchange.com/q/3179834 and links therein 
for some hints.

Duncan Murdoch

Thanks, I'll check it out. I ran the simulation and I got:
```

t = 1, N = 20,000
t = 2, N = 40,000
t = 3, N = 80,000
t = 4, N = 160,000
t = 5, N = 320,000
t = 6, N = 640,000
t = 7, N = 1,280,000
```
Hence the answer is t=6.{...} but the problem is to get that
fractional value. Would be possible to use some kind of interpolation?
I have the known Xs (the t values), the known Ys (the Nt), I need to
get x when y is 10?

On Sun, Jan 24, 2021 at 9:40 PM Duncan Murdoch <murdoch.duncan at
gmail.com> wrote:
>
> On 24/01/2021 2:57 p.m., Luigi Marongiu 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?
>
> There are two answers:
>
> The brute force answer is just to try it and count how far you need to
> go.  This is really simple, but really inefficient.
>
> The faster and more elegant way is to solve the recursive relation for
> an explicit solution.  You've got a quadratic recurrence relation;
> there's no general solution to those, but there are solutions in
special
> cases.  See https://math.stackexchange.com/q/3179834 and links therein
> for some hints.
>
> Duncan Murdoch


-- 
Best regards,
Luigi