fork(1) download 
  1. # your code goes here
  2. likelihood=function(x,y){
  3.   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))*
  4.     (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)*
  5.     (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)*
  6.     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))*
  7.     (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)*
  8.     (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))
  9. }
  10.  
  11. mu=seq(10.1,10.2,0.001)
  12. sigma=seq(0.5,0.55,0.001)
  13. z=outer(mu,sigma,likelihood)
  14. ans=contour(mu,sigma,z,label=seq(0,1,by=.1))
  15.  
  16. #求mle
  17. optimize(likelihood,c(0,20),maximum=T)
  18. nlminb(start = 0, obj = likelihood , lower = 0, upper = 1000)$par
  19.  
  20. #微分
  21. likelihood2=expression(
  22.   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))*
  23.     (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)*
  24.     (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)*
  25.     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))*
  26.     (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)*
  27.     (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))
  28. )
  29. deriv(likelihood2, namevec = "x")
  30.  
  31.  
  32. relative=function(x,y){
  33.   likelihood(x,0.53)/likelihood(10.14,0.53) 
  34. }
  35. x=seq(9,12,0.001)
  36. plot(x,relative(x),type="l")
  37. abline(h=0.147,col=2)
  38.  
  39.  
  40. relative2=function(x,y){
  41.   f=likelihood(x,0.5301)/likelihood(10.14,0.5301) 
  42.   return(f-0.147)
  43. }
  44. uniroot(relative2,c(9.5,10))$root
  45. uniroot(relative2,c(10.3,10.5))$root
  46.  
  47. abline(v=9.9253593,col=2)
  48. abline(v=10.37615,col=2)
  49.  
Success #stdin #stdout #stderr 0.26s 43892KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
stderr 
Error in dlnorm(6700, x, y) : argument "y" is missing, with no default
Calls: optimize -> <Anonymous> -> f -> dlnorm
Execution halted
https://ideone.com/V6jzsf
language:
R (R 3.5.2)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!