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

Luigi Marongiu m@rong|u@|u|g| @end|ng |rom gm@||@com
Mon Jan 25 15:20:50 CET 2021


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 using 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



More information about the R-help mailing list