function [logdensity, score, prediction, second_raw_moment] = Poisson(theta, dat
a)
	y = nth(data,1);
	x = nth(data,2);
	k = columns(x);
	lambda = exp(x*theta);
	logdensity = -lambda  +  y .* (x*theta) - mc_lgamma(y+1);
	score = y - lambda;
	score = dmult(score,x);
	prediction = lambda;
	second_raw_moment = prediction + prediction .^2;

endfunction