--- a/inst/poisson_moments.m
+++ b/inst/poisson_moments.m
@@ -22,5 +22,5 @@
 	w = data(:, k+2:columns(data));
 	lambda = exp(x*theta);
 	e = y ./ lambda - 1;
-	m = dmult(e, w);
+	m = diag(e) * w;
 endfunction