[R] Need help

J C Nash pro|jcn@@h @end|ng |rom gm@||@com
Thu Aug 4 17:32:07 CEST 2022 

FWIW, package Rmpfr has unirootR which Martin Maechler included after I coded a version of zeroin
in pure R at UseR2011. (There was a long session in biostats that was a bit outside my area, and ...)

Sometimes useful when high precision necessary.

JN


On 2022-08-04 11:24, Rui Barradas wrote:
> Hello,
> 
> Like others have said, the root does depend on the seed but not by much.
> Here is a loop finding the roots for seeds from 1 to 10,000.
> 
> I have defined a variable n == 1000 instead of having the calls to rnorm depend on a hard-coded constant.
> 
> 
> 
> ff <- function(zz){
>    inner = vector("numeric", length = s)
>    for(k in 1:s){
>      inner[k]=(1- lam*((1+b[k]*((zz-thr)/a[k]))^(-1/b[k])))
>    }
>    answer = mean(inner) - (1 - (1/r))
>    return(answer)
> }
> 
> n <- 1000
> s <- n
> lam <- 0.15
> thr <- 70
> r <- 10
> 
> up_lims <- 112.5   # found by trial and error to be nice most
>                     # of the time, see mean(ok) below
> out_list <- vector("list", length = 1e4)
> 
> for(i in seq_along(out_list)) {
>    # give the user feedback
>    if(i %% 100 == 0) print(i)
> 
>    set.seed(i)
>    a <- rnorm(n, 110, 5)
>    b <- rnorm(n, -0.3, 0.4)
> 
>    out_list[[i]] <- tryCatch(
>      uniroot(ff, lower = 0, upper = up_lims),
>      error = function(e) e
>    )
> }
> 
> ok <- !sapply(out_list, inherits, 'error')
> mean(ok)
> #[1] 0.9978
> 
> root <- sapply(out_list[ok], '[[', 'root')
> hist(root)
> 
> 
> 
> Hope this helps,
> 
> Rui Barradas
> 
> Às 15:37 de 04/08/2022, Ebert,Timothy Aaron escreveu:
>> And one last edit.
>> Rstudio Environment shows out$root as 112. However, if I then print out$root I get 111.5303 which is much closer to 
>> the right answer.  Entering values for ff, I can work out the answer to slightly more than 111.53030196.
>>
>> Regards,
>> Tim
>>
>> -----Original Message-----
>> From: Ebert,Timothy Aaron <tebert using ufl.edu>
>> Sent: Thursday, August 4, 2022 10:28 AM
>> To: Ebert,Timothy Aaron <tebert using ufl.edu>; Md. Moyazzem Hossain <hossainmm using juniv.edu>; J C Nash <profjcnash using gmail.com>
>> Cc: r-help mailing list <r-help using r-project.org>
>> Subject: RE: [R] Need help
>>
>> I wish I could edit my earlier reply.
>> ff(144.43) returns 0.03142121, but ff(144.44) returns NaN.
>> ff(12) returns -0.1468287
>>
>> out=uniroot(ff, lower=12, upper =144) does not have errors.
>> out is a list of 5.
>> $root is 112
>> $f.root is 1.09e-08
>> $iter is 5
>> $Init.it is NA
>> $estim.prec is 6.1e-05
>>
>> Is this working?
>>
>> Regards,
>> Tim
>>
>> -----Original Message-----
>> From: R-help <r-help-bounces using r-project.org> On Behalf Of Ebert,Timothy Aaron
>> Sent: Thursday, August 4, 2022 10:09 AM
>> To: Md. Moyazzem Hossain <hossainmm using juniv.edu>; J C Nash <profjcnash using gmail.com>
>> Cc: r-help mailing list <r-help using r-project.org>
>> Subject: Re: [R] Need help
>>
>> [External Email]
>>
>> A few issues
>> 1) What does this function do? If you describe the problem and goal there might be a better answer. We appreciate 
>> seeing that you have attempted a solution, but you have trapped us all in blindly following your attempt.
>>
>> 2) This is an error checking phase of programming. Setting the seed is a good idea. Just remember to unset the seed 
>> after finishing the debugging phase.
>>
>> 3) In uniroot you call the function, but the function expects an argument (zz) that is not provided, and there is no 
>> default.
>>
>> 4) I set.seed(42), and entered ff(12) getting an answer of -0.1468287. Is this the expected result?
>>
>> Regards,
>> Tim
>>
>> -----Original Message-----
>> From: R-help <r-help-bounces using r-project.org> On Behalf Of Md. Moyazzem Hossain
>> Sent: Thursday, August 4, 2022 9:50 AM
>> To: J C Nash <profjcnash using gmail.com>
>> Cc: r-help mailing list <r-help using r-project.org>
>> Subject: Re: [R] Need help
>>
>> [External Email]
>>
>> Dear JN,
>>
>> Thanks.
>>
>> I do not check whether the function actually crosses zero or not. However, by assumption, the value would be greater 
>> than zero.
>>
>> Hossain
>>
>> On Thu, Aug 4, 2022 at 2:40 PM J C Nash <profjcnash using gmail.com> wrote:
>>
>>> Have you checked that your function actually crosses zero?
>>>
>>> You should also set a seed if you want a reproducible result.
>>>
>>> JN
>>>
>>> On 2022-08-04 09:30, Md. Moyazzem Hossain wrote:
>>>> Dear R Experts,
>>>>
>>>> I hope that you are doing well.
>>>>
>>>> I am facing a problem to find out the value of the following
>>>> function. I need help in this regard.
>>>>
>>>> #####
>>>> a=rnorm(1000, 110, 5)
>>>> b = rnorm(1000, -0.3, 0.4)
>>>> s = length(a)
>>>> lam=0.15
>>>> thr=70
>>>> r= 10
>>>>
>>>> ff = function(zz){
>>>>     inner = vector("numeric", length = s)
>>>>        for(k in 1:s){
>>>>         inner[k]=(1- lam*((1+b[k]*((zz-thr)/a[k]))^(-1/b[k])))
>>>>             }
>>>>     answer = mean(inner)- (1- (1/r))
>>>>     return(answer)
>>>>     }
>>>> ########
>>>> out=uniroot(ff, lower = 0, upper = 10000 )$root out
>>>>
>>>> ########### Error ########
>>>> Error in uniroot(ff, lower = 0, upper = 10000) :
>>>>     f.upper = f(upper) is NA
>>>>
>>>> Please help me. Thanks in advance.
>>>>
>>>> Take care.
>>>>
>>>> Hossain
>>>>
>>>
>>
>>
>> -- 
>> Best Regards,
>> Md. Moyazzem Hossain
>> Associate Professor
>> Department of Statistics
>> Jahangirnagar University
>> Savar, Dhaka-1342, Bangladesh
>> Website: 
>> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.juniv.edu%2Fteachers%2Fhossainmm&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=3GVM6FLCtmP9oKjjtVzE0LZpI6LL8ldXknodpFQJUbM%3D&reserved=0 
>>
>> Research: *[image: Google Scholar]
>> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fscholar.google.com%2Fcitations%3Fhl%3Den%26user%3D-U03XCgAAAAJ&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=8EfuEse2%2FnvwuZqvKskGXN4tLnM%2FDTA8XsU7nMma04c%3D&reserved=0>* 
>> | *ResearchGate
>> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.researchgate.net%2Fprofile%2FMd_Hossain107&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=IRK69rC%2F251zGT0Kpdo8jfHI1gkQhwoHxzYmozwiYKM%3D&reserved=0>* 
>> | *ORCID iD
>> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Forcid.org%2F0000-0003-3593-6936&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=XzQ51GyGR1tEB2iOG90cX9Np8pyeXSDBeeUpC0wDfEc%3D&reserved=0>* 
>>
>>
>>          [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-help&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=VD0EWkp1hKz%2FEKuJsGAsegL9K33y0S1L2O2jDzM6EsM%3D&reserved=0 
>>
>> PLEASE do read the posting guide 
>> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r-project.org%2Fposting-guide.html&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=9wjGhM9eZTIbW6ANRoxzvEI%2BH6tHeSLbGwi7fuACKmE%3D&reserved=0 
>>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>> ______________________________________________
>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-help&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=VD0EWkp1hKz%2FEKuJsGAsegL9K33y0S1L2O2jDzM6EsM%3D&reserved=0 
>>
>> PLEASE do read the posting guide 
>> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r-project.org%2Fposting-guide.html&data=05%7C01%7Ctebert%40ufl.edu%7C92190c03aa824920e8a108da76258618%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637952200791209718%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=9wjGhM9eZTIbW6ANRoxzvEI%2BH6tHeSLbGwi7fuACKmE%3D&reserved=0 
>>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>> ______________________________________________
>> 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