What you get and what you should be getting: checking numerical code

[Update 24.08.12: added missing alias following comment by brobar]

Whenever I write numerical code I spend half my time debugging my algebra, painstakingly uncovering one sign mistake after another in my calculations. Usually I have computed by hand the gradient or the integral of some nasty function, and I have to check it against a numerical estimate (for example returned by the wonderful numDeriv package). I have written a small utility called cmpplot, for “comparison plot”, that allows me to quickly plot what I’m getting from my hand-computed function against what I should be getting (the numerical estimate). It’s rather useful if you need to diagnose sign mistakes, forgotten factors, etc. It’s too small to be a R package so I’m just going to put the code below, with an example.

Let’s say that after lengthy calculations we have decided that the derivative of cos(x) is sin(x):

fun <- function(x) cos(x)
dfun <- function(x) sin(x)

Although quite confident in our result, we decide to check against a numerical derivative anyway:

The points should all fall on the diagonal line, and most obviously they don’t. cmpplot also displays the result of a regression of d.th on d.num, which tells us we’ve made a sign mistake and that the derivative of cos(x) is actually -sin(x).

in cases like this calling cmpplot with argument “type = 2″ provides a better way of visualising what’s going on:

cmpplot(d.num,d.th,type=2)

The black dots should all align with the blue dots and only half of them do so.

cmpplot also has methods for matrices: this time we take cos(x) as a function that takes a vector argument and returns a vector argument, and we want to compute its Jacobian. Again we have made a mistake in our theoretical calculations:

Related

4 comments

i have a question.
when i entered ” cmpplot(d.num,d.th) “, i found this message:
‘ error: function “cmpplot” is not found ‘
??cmpplot and help(cmpplot) also do not help for me.
could i know what i have to do?