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 at ufl.edu>
> Sent: Thursday, August 4, 2022 10:28 AM
> To: Ebert,Timothy Aaron <tebert at ufl.edu>; Md. Moyazzem Hossain
<hossainmm at juniv.edu>; J C Nash <profjcnash at gmail.com>
> Cc: r-help mailing list <r-help at 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 at r-project.org> On Behalf Of
Ebert,Timothy Aaron
> Sent: Thursday, August 4, 2022 10:09 AM
> To: Md. Moyazzem Hossain <hossainmm at juniv.edu>; J C Nash
<profjcnash at gmail.com>
> Cc: r-help mailing list <r-help at 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 at r-project.org> On Behalf Of Md.
Moyazzem Hossain
> Sent: Thursday, August 4, 2022 9:50 AM
> To: J C Nash <profjcnash at gmail.com>
> Cc: r-help mailing list <r-help at 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 at 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 at 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 at 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 at 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.