{"action":"opened","issue":{"url":"https://api.github.com/repos/kaskr/adcomp/issues/274","repository_url":"https://api.github.com/repos/kaskr/adcomp","labels_url":"https://api.github.com/repos/kaskr/adcomp/issues/274/labels{/name}","comments_url":"https://api.github.com/repos/kaskr/adcomp/issues/274/comments","events_url":"https://api.github.com/repos/kaskr/adcomp/issues/274/events","html_url":"https://github.com/kaskr/adcomp/issues/274","id":331028587,"node_id":"MDU6SXNzdWUzMzEwMjg1ODc=","number":274,"title":"Why autodiff::hessian always gives wrong results?","user":{"login":"Zednathan","id":38760310,"node_id":"MDQ6VXNlcjM4NzYwMzEw","avatar_url":"https://avatars2.githubusercontent.com/u/38760310?v=4","gravatar_id":"","url":"https://api.github.com/users/Zednathan","html_url":"https://github.com/Zednathan","followers_url":"https://api.github.com/users/Zednathan/followers","following_url":"https://api.github.com/users/Zednathan/following{/other_user}","gists_url":"https://api.github.com/users/Zednathan/gists{/gist_id}","starred_url":"https://api.github.com/users/Zednathan/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/Zednathan/subscriptions","organizations_url":"https://api.github.com/users/Zednathan/orgs","repos_url":"https://api.github.com/users/Zednathan/repos","events_url":"https://api.github.com/users/Zednathan/events{/privacy}","received_events_url":"https://api.github.com/users/Zednathan/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":0,"created_at":"2018-06-11T02:40:29Z","updated_at":"2018-06-11T02:40:29Z","closed_at":null,"author_association":"NONE","body":"I wrote a code, but autodiff::hessian() always gives wrong results when I change the parameters values, but all the other things are correct. Could you please tell me whether I used autodiff::hessian() correctly? I made it bold text in the code. Thank you very much.\r\nMy code is as follows:\r\n#include <TMB.hpp>\r\n# include <iostream >\r\n/* The following is the culculation for a single peak */\r\ntemplate<class Type>\r\nstruct joint_nll {\r\n  /* Data and parameter objects for a single peak: */\r\n  Type log_length0;\r\n  Type log_age0;\r\n  Type log_mu0;\r\n  Type sig_log_len_me;\r\n  Type sig_log_age_me;\r\n  Type mix_mu0;\r\n  Type mix_sigma0;\r\n  /* Constructor */\r\n  joint_nll(Type log_length0_, Type log_age0_, Type log_mu0_, Type sig_log_len_me_,\r\n            Type sig_log_age_me_, Type mix_mu0_, Type mix_sigma0_) :\r\n    log_length0(log_length0_), log_age0(log_age0_), log_mu0(log_mu0_), sig_log_len_me(sig_log_len_me_),\r\n    sig_log_age_me(sig_log_age_me_), mix_mu0(mix_mu0_), mix_sigma0(mix_sigma0_) {}\r\n  /* Evaluate the negative joint log-likelihood as function of the\r\n     random effects */\r\n  template <typename T>\r\n  T operator()(vector<T> u) {\r\n    T u0 = u(0);\r\n//    T res = -T(dnorm(log_length0,log_mu0,sig_log_len_me,true));\r\n    T res = -dnorm(T(log_length0),T(log_mu0),T(sig_log_len_me),true);\r\n    res -= dnorm(T(log_age0),u0,T(sig_log_age_me),true);\r\n    res -= dnorm(u0,T(mix_mu0),T(mix_sigma0),true);\r\n    return res;\r\n  }\r\n};\r\n\r\ntemplate<class Type>\r\nType objective_function<Type>::operator() ()\r\n{\r\nDATA_VECTOR(age);\r\nDATA_VECTOR(length);\r\nDATA_VECTOR(log_age);\r\nDATA_VECTOR(log_length);\r\nDATA_SCALAR(sig_log_len_me);\r\nDATA_SCALAR(sig_log_age_me);\r\n\r\nPARAMETER(log_Linf);\r\nPARAMETER(log_k);\r\nPARAMETER(ao);\r\nPARAMETER_VECTOR(lambda);\r\nPARAMETER_VECTOR(mix_mu_parm);\r\nPARAMETER_VECTOR(log_mix_sigma);\r\nPARAMETER_VECTOR(log_true_age);\r\n\r\nint n = age.size();\r\nint nmix = lambda.size();\r\nType Linf = exp(log_Linf);\r\nType k = exp(log_k);\r\nvector <Type> mix_sigma = exp(log_mix_sigma);\r\nvector <Type> true_age = exp(log_true_age);\r\nvector <Type> mu = Linf*(Type(1.0) - exp(-k*(true_age-ao)));\r\nvector <Type> log_mu = log(mu);\r\nvector <Type> p(nmix+1);\r\nvector <Type> mix_mu(nmix+1);\r\nvector <Type> log_p(nmix+1);\r\nType lterm = 0.0;\r\nfor(int i = 0;i <nmix;++i){\r\nlterm += log(1 + exp(lambda(i)));\r\nlog_p(i) = lambda(i) - lterm;}\r\nlog_p(nmix) = -lterm;\r\np = exp(log_p);\r\nmix_mu(0) = mix_mu_parm(0);\r\nfor(int i = 1;i <= nmix;++i){\r\nmix_mu(i) = mix_mu(i-1) + exp(mix_mu_parm(i));\r\n}\r\n\r\nType nll = 0.0;\r\nfor(int j = 0;j<=nmix;++j){\r\nfor(int i = 0;i<n;++i){\r\nnll -= dnorm(log_length(i),log_mu(j*n+i),sig_log_len_me,true);\r\nnll -= dnorm(log_age(i),log_true_age(j*n+i),sig_log_age_me,true);\r\nnll -= dnorm(log_true_age(j*n+i),mix_mu(j),mix_sigma(j),true);\r\n}\r\n}\r\n\r\nType nll_true = 0.0;\r\nfor(int i = 0;i<n;++i){\r\nType nlli = 0.0;\r\nfor(int j = 0;j<=nmix;++j){\r\n  // Construct joint negative log-likelihood\r\n  joint_nll<Type> jnll(log_length(i),log_age(i),log_mu(j*n+i),sig_log_len_me,sig_log_age_me,mix_mu(j),mix_sigma(j));\r\n  // Calculate Laplace approx\r\n  vector<Type> u(2);\r\n  u(0) = log_true_age(j*n+i);\r\n**// here the autodiff::hessian does not give correct results**\r\n  matrix<Type> H = autodiff::hessian(jnll, u);\r\n  Type res = .5 * log(H(0,0)) + jnll(u);\r\n  res -= .5 * log(2.0 * M_PI);\r\n  nlli += p(j)*exp(-res);\r\n}\r\nnll_true -= log(nlli);\r\n}\r\n\r\nREPORT(nll_true);\r\nREPORT(Linf);  \r\nREPORT(k);    \r\nREPORT(ao);   \r\n \r\nreturn nll;\r\n}\r\n"}}
{"actor":{"display_login":"Zednathan"}}