How to find when a value is reached given a function?

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

How to find when a value is reached given a function?

Luigi
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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: How to find when a value is reached given a function?

Duncan Murdoch-2
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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: How to find when a value is reached given a function?

Luigi
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 <[hidden email]> 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: How to find when a value is reached given a function?

Luigi
If I run this:
```
Y = c(10000, 20000, 40000, 80000, 160000, 320000, 640000, 1280000)
X = 0:7
plot(Y~X, log='y')
model <- lm(log10(Y) ~ X)
abline(model)
predict(model, data.frame(Y=log10(1000000)))
```
I get a funny answer:
```
      1       2       3       4       5       6       7       8
4.00000 4.30103 4.60206 4.90309 5.20412 5.50515 5.80618 6.10721
Warning message:
'newdata' had 1 row but variables found have 8 rows
```
but:
```
> data.frame(Y=log10(1000000))
  Y
1 6
```
What is the correct use?
Thank you

On Mon, Jan 25, 2021 at 3:20 PM Luigi Marongiu <[hidden email]> wrote:

>
> 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 <[hidden email]> 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



--
Best regards,
Luigi

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: How to find when a value is reached given a function?

aBBy Spurdle, ⍺XY
In reply to this post by Luigi
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 <[hidden email]> 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
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.