Re: [R] Query on finding root

[Ben Bolker]

(I mean pdavies)

On Mon, Aug 28, 2023, 7:52 AM Ben Bolker  wrote:

> I would probably use the built in qdavies() function...
>
> On Mon, Aug 28, 2023, 7:48 AM Leonard Mada via R-help <
> r-help@r-project.org> wrote:
>
>> Dear R-Users,
>>
>> Just out of curiosity:
>> Which of the 2 methods is the better one?
>>
>> The results seem to differ slightly.
>>
>>
>> fun = function(u){((26104.50*u^0.03399381)/((1-u)^0.107)) - 28353.7}
>>
>> uniroot(fun, c(0,1))
>> # 0.6048184
>>
>> curve(fun(x), 0, 1)
>> abline(v=0.3952365, col="red")
>> abline(v=0.6048184, col="red")
>> abline(h=0, col="blue")
>>
>>
>>
>> fun = function(u){ (0.03399381*log(u) - 0.107*log(1-u)) -
>> log(28353.7/26104.50) }
>> fun = function(u){ (0.03399381*log(u) - 0.107*log1p(-u)) -
>> log(28353.7/26104.50) }
>>
>> uniroot(fun, c(0,1))
>> # 0.6047968
>>
>> curve(fun(x), 0, 1)
>> abline(v=0.3952365, col="red")
>> abline(v=0.6047968, col="red")
>> abline(h=0, col="blue")
>>
>> Sincerely,
>>
>> Leonard
>>
>> __
>> R-help@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.
>>
>

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Query on finding root

[Ben Bolker]

I would probably use the built in qdavies() function...

On Mon, Aug 28, 2023, 7:48 AM Leonard Mada via R-help 
wrote:

> Dear R-Users,
>
> Just out of curiosity:
> Which of the 2 methods is the better one?
>
> The results seem to differ slightly.
>
>
> fun = function(u){((26104.50*u^0.03399381)/((1-u)^0.107)) - 28353.7}
>
> uniroot(fun, c(0,1))
> # 0.6048184
>
> curve(fun(x), 0, 1)
> abline(v=0.3952365, col="red")
> abline(v=0.6048184, col="red")
> abline(h=0, col="blue")
>
>
>
> fun = function(u){ (0.03399381*log(u) - 0.107*log(1-u)) -
> log(28353.7/26104.50) }
> fun = function(u){ (0.03399381*log(u) - 0.107*log1p(-u)) -
> log(28353.7/26104.50) }
>
> uniroot(fun, c(0,1))
> # 0.6047968
>
> curve(fun(x), 0, 1)
> abline(v=0.3952365, col="red")
> abline(v=0.6047968, col="red")
> abline(h=0, col="blue")
>
> Sincerely,
>
> Leonard
>
> __
> R-help@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.
>

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Query on finding root

[Leonard Mada via R-help]

Dear R-Users,

Just out of curiosity:
Which of the 2 methods is the better one?

The results seem to differ slightly.


fun = function(u){((26104.50*u^0.03399381)/((1-u)^0.107)) - 28353.7}

uniroot(fun, c(0,1))
# 0.6048184

curve(fun(x), 0, 1)
abline(v=0.3952365, col="red")
abline(v=0.6048184, col="red")
abline(h=0, col="blue")



fun = function(u){ (0.03399381*log(u) - 0.107*log(1-u)) - 
log(28353.7/26104.50) }
fun = function(u){ (0.03399381*log(u) - 0.107*log1p(-u)) - 
log(28353.7/26104.50) }


uniroot(fun, c(0,1))
# 0.6047968

curve(fun(x), 0, 1)
abline(v=0.3952365, col="red")
abline(v=0.6047968, col="red")
abline(h=0, col="blue")

Sincerely,

Leonard

__
R-help@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.

Re: [R] Query on finding root

[Ben Bolker]

   This doesn't look like homework to me -- too specific.  The posting 
guide  says that the list 
is not intended for "Basic statistics and classroom homework" -- again, 
this doesn't seem to fall into that category.


  tl;dr, I think the difference between the two approaches is just 
whether the lower or upper tail is considered (i.e., adding lower.tail = 
FALSE to pdavies(), or more simply taking (1-x), makes the two answers 
agree).


  If you look at the source code for pdavies() you'll see that it's 
essentially doing the same uniroot() calculation that you are.


## Q(u)=(c*u^lamda1)/((1-u)^lamda2)
mean <- 28353.7 # mean calculated from data
lambda1 <- .03399381 # estimates c, lambda1 and lambda2 calculated from data
lambda2 <- .107
c <- 26104.50
library(Davies)# using package
params <- c(c,lambda1,lambda2)
u <- pdavies(x = mean, params = params, lower.tail = FALSE)
u
fun <- function(u) {
with(as.list(params), (c*u^lambda1)/((1-u)^lambda2)) - mean
}
curve(fun, from = 0.01, to = 1)
uniroot <- uniroot(fun,c(0.01,1))
abline(h = 0)
uniroot$root



On 2023-08-27 5:40 p.m., Rolf Turner wrote:


On Fri, 25 Aug 2023 22:17:05 +0530
ASHLIN VARKEY  wrote:


Sir,


Please note that r-help is a mailing list, not a knight! ️


I want to solve the equation Q(u)=mean, where Q(u) represents the
quantile function. Here my Q(u)=(c*u^lamda1)/((1-u)^lamda2), which is
the quantile function of Davies (Power-pareto) distribution.  Hence I
want to solve , *(c*u^lamda1)/((1-u)^lamda2)=28353.7(Eq.1)*
where lamda1=0.03399381, lamda2=0.107 and c=26104.50. When I used
the package 'Davies' and solved Eq 1, I got the answer u=0.3952365.
But when I use the function  'uniroot' to solve the Eq.1, I got a
different answer which is  u=0.6048157.  Why did this difference
happen?  Which is the correct method to solve Eq.1. Using the value
of *u *from the first method my further calculation was nearer to
empirical values.  The R-code I used is herewith. Kindly help me to
solve this issue.

R-code
Q(u)=(c*u^lamda1)/((1-u)^lamda2)
mean=28353.7 # mean calculated from data
lamda1=.03399381 # estimates c, lamda1 and lamda2 calculated from data
lamda2=.107
c=26104.50
library(Davies)# using package
params=c(c,lamda1,lamda2)
u=pdavies(28353.7,params)
u
fun=function(u){((26104.50*u^0.03399381)/((1-u)^0.107))-28353.7}
uniroot= uniroot(fun,c(0.01,1))
uniroot


As Prof. Nash has pointed out, this looks like homework.

Some general advice:  graphics can be very revealing, and are easy to
effect in R.  Relevant method: plot.function(); relevant utility:
abline().  Look at the help for these.

cheers,

Rolf Turner



--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
(Acting) Graduate chair, Mathematics & Statistics
> E-mail is sent at my convenience; I don't expect replies outside of 
working hours.


__
R-help@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.

Re: [R] Query on finding root

[Rolf Turner]

On Fri, 25 Aug 2023 22:17:05 +0530
ASHLIN VARKEY  wrote:

> Sir,

Please note that r-help is a mailing list, not a knight! ️

> I want to solve the equation Q(u)=mean, where Q(u) represents the
> quantile function. Here my Q(u)=(c*u^lamda1)/((1-u)^lamda2), which is
> the quantile function of Davies (Power-pareto) distribution.  Hence I
> want to solve , *(c*u^lamda1)/((1-u)^lamda2)=28353.7(Eq.1)*
> where lamda1=0.03399381, lamda2=0.107 and c=26104.50. When I used
> the package 'Davies' and solved Eq 1, I got the answer u=0.3952365.
> But when I use the function  'uniroot' to solve the Eq.1, I got a
> different answer which is  u=0.6048157.  Why did this difference
> happen?  Which is the correct method to solve Eq.1. Using the value
> of *u *from the first method my further calculation was nearer to
> empirical values.  The R-code I used is herewith. Kindly help me to
> solve this issue.
> 
> R-code
> Q(u)=(c*u^lamda1)/((1-u)^lamda2)
> mean=28353.7 # mean calculated from data
> lamda1=.03399381 # estimates c, lamda1 and lamda2 calculated from data
> lamda2=.107
> c=26104.50
> library(Davies)# using package
> params=c(c,lamda1,lamda2)
> u=pdavies(28353.7,params)
> u
> fun=function(u){((26104.50*u^0.03399381)/((1-u)^0.107))-28353.7}
> uniroot= uniroot(fun,c(0.01,1))
> uniroot

As Prof. Nash has pointed out, this looks like homework.

Some general advice:  graphics can be very revealing, and are easy to
effect in R.  Relevant method: plot.function(); relevant utility:
abline().  Look at the help for these.

cheers,

Rolf Turner

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Stats. Dep't. (secretaries) phone:
 +64-9-373-7599 ext. 89622
Home phone: +64-9-480-4619

__
R-help@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.

Re: [R] Query on finding root

[J C Nash]

Homework?

On 2023-08-25 12:47, ASHLIN VARKEY wrote:

Sir,
I want to solve the equation Q(u)=mean, where Q(u) represents the quantile
function. Here my Q(u)=(c*u^lamda1)/((1-u)^lamda2), which is the quantile
function of Davies (Power-pareto) distribution.  Hence I want to solve ,
*(c*u^lamda1)/((1-u)^lamda2)=28353.7(Eq.1)*
where lamda1=0.03399381, lamda2=0.107 and c=26104.50. When I used the
package 'Davies' and solved Eq 1, I got the answer u=0.3952365. But when I
use the function  'uniroot' to solve the Eq.1, I got a different answer
which is  u=0.6048157.  Why did this difference happen?  Which is the
correct method to solve Eq.1. Using the value of *u *from the first method
my further calculation was nearer to empirical values.  The R-code I used
is herewith. Kindly help me to solve this issue.

R-code
Q(u)=(c*u^lamda1)/((1-u)^lamda2)
mean=28353.7 # mean calculated from data
lamda1=.03399381 # estimates c, lamda1 and lamda2 calculated from data
lamda2=.107
c=26104.50
library(Davies)# using package
params=c(c,lamda1,lamda2)
u=pdavies(28353.7,params)
u
fun=function(u){((26104.50*u^0.03399381)/((1-u)^0.107))-28353.7}
uniroot= uniroot(fun,c(0.01,1))
uniroot

[[alternative HTML version deleted]]

__
R-help@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.


__
R-help@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.

[R] Query on finding root

[ASHLIN VARKEY]

Sir,
I want to solve the equation Q(u)=mean, where Q(u) represents the quantile
function. Here my Q(u)=(c*u^lamda1)/((1-u)^lamda2), which is the quantile
function of Davies (Power-pareto) distribution.  Hence I want to solve ,
*(c*u^lamda1)/((1-u)^lamda2)=28353.7(Eq.1)*
where lamda1=0.03399381, lamda2=0.107 and c=26104.50. When I used the
package 'Davies' and solved Eq 1, I got the answer u=0.3952365. But when I
use the function  'uniroot' to solve the Eq.1, I got a different answer
which is  u=0.6048157.  Why did this difference happen?  Which is the
correct method to solve Eq.1. Using the value of *u *from the first method
my further calculation was nearer to empirical values.  The R-code I used
is herewith. Kindly help me to solve this issue.

R-code
Q(u)=(c*u^lamda1)/((1-u)^lamda2)
mean=28353.7 # mean calculated from data
lamda1=.03399381 # estimates c, lamda1 and lamda2 calculated from data
lamda2=.107
c=26104.50
library(Davies)# using package
params=c(c,lamda1,lamda2)
u=pdavies(28353.7,params)
u
fun=function(u){((26104.50*u^0.03399381)/((1-u)^0.107))-28353.7}
uniroot= uniroot(fun,c(0.01,1))
uniroot

[[alternative HTML version deleted]]

__
R-help@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.