diff --git a/matlab/lpdfgam.m b/matlab/lpdfgam.m
index 0a8505960..1246f2a9e 100644
--- a/matlab/lpdfgam.m
+++ b/matlab/lpdfgam.m
@@ -1,4 +1,4 @@
-function  ldens = lpdfgam(x,a,b);
+function  [ldens,Dldens,D2ldens] = lpdfgam(x,a,b);
 % Evaluates the logged GAMMA PDF at x.
 %
 % INPUTS     
@@ -37,4 +37,22 @@ if length(a)==1
     ldens(idx) = -gammaln(a) - a*log(b) + (a-1)*log(x(idx)) - x(idx)/b ;
 else
     ldens(idx) = -gammaln(a(idx)) - a(idx).*log(b(idx)) + (a(idx)-1).*log(x(idx)) - x(idx)./b(idx) ;
+end
+
+
+
+if nargout >1 
+    if length(a)==1
+        Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ;
+    else
+        Dldens(idx) = (a(idx)-1)./(x(idx)) - ones(length(idx),1)./b(idx) ;
+    end
+end
+
+if nargout == 3
+    if length(a)==1
+        D2ldens(idx) = -(a-1)./(x(idx)).^2;
+    else
+        D2ldens(idx) = -(a(idx)-1)./(x(idx)).^2;
+    end
 end
\ No newline at end of file
diff --git a/matlab/lpdfgbeta.m b/matlab/lpdfgbeta.m
index 33ff3df48..1906161e5 100644
--- a/matlab/lpdfgbeta.m
+++ b/matlab/lpdfgbeta.m
@@ -1,4 +1,4 @@
-function ldens = lpdfgbeta(x,a,b,aa,bb);
+function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb);
 % Evaluates the logged BETA PDF at x. 
 %
 % INPUTS 
@@ -38,4 +38,22 @@ if length(a)==1
     ldens(idx) = -betaln(a,b) + (a-1)*log(x(idx)-aa) + (b-1)*log(bb-x(idx)) - (a+b-1)*log(bb-aa) ;
 else
     ldens(idx) = -betaln(a(idx),b(idx)) + (a(idx)-1).*log(x(idx)-aa(idx)) + (b(idx)-1).*log(bb(idx)-x(idx)) - (a(idx)+b(idx)-1).*log(bb(idx)-aa(idx));
+end
+
+
+if nargout >1 
+    if length(a)==1
+        Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ;
+    else
+        Dldens(idx) = (a(idx)-1)./(x(idx)-aa(idx)) - (b(idx)-1)./(bb(idx)-x(idx));
+    end
+end
+
+
+if nargout == 3 
+    if length(a)==1
+        D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 + (b-1)./(bb-x(idx)).^2 ;
+    else
+        D2ldens(idx) = -(a(idx)-1)./(x(idx)-aa(idx)).^2 + (b(idx)-1)./(bb(idx)-x(idx)).^2;
+    end
 end
\ No newline at end of file
diff --git a/matlab/lpdfig1.m b/matlab/lpdfig1.m
index 9c17277a8..bb0b54078 100644
--- a/matlab/lpdfig1.m
+++ b/matlab/lpdfig1.m
@@ -1,4 +1,4 @@
-function ldens = lpdfig1(x,s,nu)
+function [ldens,Dldens,D2ldens] = lpdfig1(x,s,nu)
 % Evaluates the logged INVERSE-GAMMA-1 PDF at x.
 %
 % X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution) 
@@ -41,4 +41,21 @@ if length(s)==1
     ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ;
 else
     ldens(idx) = log(2) - gammaln(.5*nu(idx)) - .5*nu(idx).*(log(2)-log(s(idx))) - (nu(idx)+1).*log(x(idx)) - .5*s(idx)./(x(idx).*x(idx)) ;
-end
\ No newline at end of file
+end
+
+if nargout >1 
+    if length(s)==1
+        Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ;
+    else
+        Dldens(idx) = - (nu(idx)+1)./(x(idx)) + s(idx)./(x(idx).^3) ;
+    end
+end
+
+if nargout == 3 
+    if length(s)==1
+        D2ldens(idx) =  (nu+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ;
+    else
+        D2ldens(idx) =  (nu(idx)+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ;
+    end
+end
+
diff --git a/matlab/lpdfig2.m b/matlab/lpdfig2.m
index 7f9558181..63e458db0 100644
--- a/matlab/lpdfig2.m
+++ b/matlab/lpdfig2.m
@@ -1,4 +1,4 @@
-function ldens = lpdfig2(x,s,nu)
+function [ldens,Dldens,D2ldens] = lpdfig2(x,s,nu)
 % Evaluates the logged INVERSE-GAMMA-2 PDF at x.
 %
 % X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) 
@@ -41,4 +41,20 @@ if length(s)==1
     ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx);
 else
     ldens(idx) = -gammaln(.5*nu(idx)) - (.5*nu(idx)).*(log(2)-log(s(idx))) - .5*(nu(idx)+2).*log(x(idx)) -.5*s(idx)./x(idx);
+end
+
+if nargout >1 
+    if length(s)==1
+        Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2;
+    else
+        Dldens(idx) = - .5*(nu(idx)+2)./(x(idx)) + .5*s(idx)./x(idx).^2;
+    end
+end
+
+if nargout == 3 
+    if length(s)==1
+        D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3;
+    else
+        D2ldens(idx) = .5*(nu(idx)+2)./(x(idx)).^2 - s(idx)./x(idx).^3;
+    end
 end
\ No newline at end of file
diff --git a/matlab/lpdfnorm.m b/matlab/lpdfnorm.m
index 0e232884a..9afdd2b6d 100644
--- a/matlab/lpdfnorm.m
+++ b/matlab/lpdfnorm.m
@@ -1,4 +1,4 @@
-function  ldens = lpdfnorm(x,a,b)
+function  [ldens,Dldens,D2ldens] = lpdfnorm(x,a,b)
 % Evaluates the logged UNIVARIATE GAUSSIAN PDF at x.
 %
 % INPUTS 
@@ -32,4 +32,12 @@ function  ldens = lpdfnorm(x,a,b)
 
 if nargin<3, b=1; end
 if nargin<2, a=0; end
-ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;
\ No newline at end of file
+ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;
+
+if nargout >1 
+    Dldens =  - (1/b)*((x-a)/b) ;
+end
+
+if nargout == 3 
+    D2ldens =  - (1/b)^2 ;
+end
\ No newline at end of file