title: calculate vector valued Hessian in R
body: <p>I want to calculate a variance-covariance matrix of parameters. The parameters are obtained by a non-linear least squares fit.</p>
  
  <pre><code>library(minpack.lm)
  library(numDeriv)
  </code></pre>
  
  <h3>variables</h3>
  
  <pre><code>t &lt;- seq(0.1,20,0.3)
  a &lt;- 20
  b &lt;- 14
  c &lt;- 0.4
  jitter &lt;- rnorm(length(t),0,0.5)
  Hobs &lt;- a+b*exp(-c*t)+jitter
  </code></pre>
  
  <h3>function def</h3>
  
  <pre><code>Hhat &lt;- function(parList, t) {parList$a + parList$b*exp(-parL
  Hhatde &lt;- function(par, t) {par[1] + par[2]*exp(-par[3]*t)}st$c*t)}
  residFun &lt;- function(par, t, observed) observed - Hhat(par,t)
  </code></pre>
  
  <h3>initial conditions</h3>
  
  <pre><code>parStart = list(a = 20, b = 10 ,c = 0.5)
  </code></pre>
  
  <h3>nls.lm</h3>
  
  <pre><code>library(minpack.lm)
  out1 &lt;- nls.lm(par = parStart, fn = residFun, observed = Hobs,
                t = t, control = nls.lm.control(nprint=0))
  </code></pre>
  
  <p>I wish to calculate manually what is given back via <code>vcov(out1)</code>
  I tried it with: but <code>sigma</code> and <code>vcov(out1)</code> which don't seem to be the same</p>
  
  <pre><code>J &lt;- jacobian(Hhatde, c(19.9508523,14.6586555,0.4066367 ), method="Richardson", 
  method.args=list(),t=t)
  sigma &lt;- solve((t(J)%*%J))
  vcov(out1)
  </code></pre>
  
  <p>now trying to do it with the hessian, I can't get it working for error message cf below</p>
  
  <h3>hessian</h3>
  
  <pre><code>H &lt;- hessian(Hhatde, x = c(19.9508523,14.6586555,0.4066367 ), method="complex", method.args=list(),t=t)
  
  Error in hessian.default(Hhatde, x = c(19.9508523, 14.6586555, 0.4066367),  : 
    Richardson method for hessian assumes a scalar valued function.
  </code></pre>
  
  <p>How do I do I get my <code>hessian()</code> to work.</p>
  
  <p>I am not very strong on the math here, hence the trial and error approach.</p>