# your code goes here
likelihood=function(x,y){
  dlnorm(6700,x,y)*(1-plnorm(6950,x,y))*(1-plnorm(7820,x,y))*(1-plnorm(8790,x,y))*dlnorm(9120,x,y)*(1-plnorm(9660,x,y))*(1-plnorm(9820,x,y))*
    (1-plnorm(11310,x,y))*(1-plnorm(11690,x,y))*(1-plnorm(11850,x,y))*(1-plnorm(11880,x,y))*(1-plnorm(12140,x,y))*dlnorm(12200,x,y)*
    (1-plnorm(12870,x,y))*dlnorm(13150,x,y)*(1-plnorm(13330,x,y))*(1-plnorm(13470,x,y))*(1-plnorm(14040,x,y))*dlnorm(14300,x,y)*
    dlnorm(17520,x,y)*(1-plnorm(17540,x,y))*(1-plnorm(17890,x,y))*(1-plnorm(18450,x,y))*(1-plnorm(18960,x,y))*(1-plnorm(18980,x,y))*
    (1-plnorm(19410,x,y))*dlnorm(20100,x,y)*(1-plnorm(20100,x,y))*(1-plnorm(20150,x,y))*(1-plnorm(20320,x,y))*dlnorm(20900,x,y)*dlnorm(22700,x,y)*
    (1-plnorm(23490,x,y))*dlnorm(26510,x,y)*(1-plnorm(27410,x,y))*dlnorm(27490,x,y)*(1-plnorm(27890,x,y))*(1-plnorm(28100,x,y))
}

mu=seq(10.1,10.2,0.001)
sigma=seq(0.5,0.55,0.001)
z=outer(mu,sigma,likelihood)
ans=contour(mu,sigma,z,label=seq(0,1,by=.1))

#求mle
optimize(likelihood,c(0,20),maximum=T)
nlminb(start = 0, obj = likelihood , lower = 0, upper = 1000)$par

#微分
likelihood2=expression(
  dlnorm(6700,x,y)*(1-plnorm(6950,x,y))*(1-plnorm(7820,x,y))*(1-plnorm(8790,x,y))*dlnorm(9120,x,y)*(1-plnorm(9660,x,y))*(1-plnorm(9820,x,y))*
    (1-plnorm(11310,x,y))*(1-plnorm(11690,x,y))*(1-plnorm(11850,x,y))*(1-plnorm(11880,x,y))*(1-plnorm(12140,x,y))*dlnorm(12200,x,y)*
    (1-plnorm(12870,x,y))*dlnorm(13150,x,y)*(1-plnorm(13330,x,y))*(1-plnorm(13470,x,y))*(1-plnorm(14040,x,y))*dlnorm(14300,x,y)*
    dlnorm(17520,x,y)*(1-plnorm(17540,x,y))*(1-plnorm(17890,x,y))*(1-plnorm(18450,x,y))*(1-plnorm(18960,x,y))*(1-plnorm(18980,x,y))*
    (1-plnorm(19410,x,y))*dlnorm(20100,x,y)*(1-plnorm(20100,x,y))*(1-plnorm(20150,x,y))*(1-plnorm(20320,x,y))*dlnorm(20900,x,y)*dlnorm(22700,x,y)*
    (1-plnorm(23490,x,y))*dlnorm(26510,x,y)*(1-plnorm(27410,x,y))*dlnorm(27490,x,y)*(1-plnorm(27890,x,y))*(1-plnorm(28100,x,y))
)
deriv(likelihood2, namevec = "x")


relative=function(x,y){
  likelihood(x,0.53)/likelihood(10.14,0.53) 
}
x=seq(9,12,0.001)
plot(x,relative(x),type="l")
abline(h=0.147,col=2)


relative2=function(x,y){
  f=likelihood(x,0.5301)/likelihood(10.14,0.5301) 
  return(f-0.147)
}
uniroot(relative2,c(9.5,10))$root
uniroot(relative2,c(10.3,10.5))$root

abline(v=9.9253593,col=2)
abline(v=10.37615,col=2)