diff options
Diffstat (limited to 'main_func/studentst.m')
-rw-r--r-- | main_func/studentst.m | 94 |
1 files changed, 47 insertions, 47 deletions
diff --git a/main_func/studentst.m b/main_func/studentst.m index 99fed1e..93e0a0a 100644 --- a/main_func/studentst.m +++ b/main_func/studentst.m @@ -1,47 +1,47 @@ -function [f,g,h,s,k] = studentst(r,k,s) -% Students T penalty with 'auto-tuning' -% -% use: -% [f,g,h,{k,{s}}] = studentst(r) - automatically fits s and k -% [f,g,h,{k,{s}}] = studentst(r,k) - automatically fits s -% [f,g,h,{k,{s}}] = studentst(r,k,s) - use given s and k -% -% input: -% r - residual as column vector -% s - scale (optional) -% k - degrees of freedom (optional) -% -% output: -% f - misfit (scalar) -% g - gradient (column vector) -% h - positive approximation of the Hessian (column vector, Hessian is a diagonal matrix) -% s,k - scale and degrees of freedom -% -% Tristan van Leeuwen, 2012. -% tleeuwen@eos.ubc.ca - -% fit both s and k -if nargin == 1 - opts = optimset('maxFunEvals',1e2); - tmp = fminsearch(@(x)st(r,x(1),x(2)),[1;2],opts); - s = tmp(1); - k = tmp(2); -end - - -if nargin == 2 - opts = optimset('maxFunEvals',1e2); - tmp = fminsearch(@(x)st(r,x,k),[1],opts); - s = tmp(1); -end - -% evaulate penalty -[f,g,h] = st(r,s,k); - - -function [f,g,h] = st(r,s,k) -n = length(r); -c = -n*(gammaln((k+1)/2) - gammaln(k/2) - .5*log(pi*s*k)); -f = c + .5*(k+1)*sum(log(1 + conj(r).*r/(s*k))); -g = (k+1)*r./(s*k + conj(r).*r); -h = (k+1)./(s*k + conj(r).*r); +function [f,g,h,s,k] = studentst(r,k,s)
+% Students T penalty with 'auto-tuning'
+%
+% use:
+% [f,g,h,{k,{s}}] = studentst(r) - automatically fits s and k
+% [f,g,h,{k,{s}}] = studentst(r,k) - automatically fits s
+% [f,g,h,{k,{s}}] = studentst(r,k,s) - use given s and k
+%
+% input:
+% r - residual as column vector
+% s - scale (optional)
+% k - degrees of freedom (optional)
+%
+% output:
+% f - misfit (scalar)
+% g - gradient (column vector)
+% h - positive approximation of the Hessian (column vector, Hessian is a diagonal matrix)
+% s,k - scale and degrees of freedom
+%
+% Tristan van Leeuwen, 2012.
+% tleeuwen@eos.ubc.ca
+
+% fit both s and k
+if nargin == 1
+ opts = optimset('maxFunEvals',1e2);
+ tmp = fminsearch(@(x)st(r,x(1),x(2)),[1;2],opts);
+ s = tmp(1);
+ k = tmp(2);
+end
+
+
+if nargin == 2
+ opts = optimset('maxFunEvals',1e2);
+ tmp = fminsearch(@(x)st(r,x,k),[1],opts);
+ s = tmp(1);
+end
+
+% evaulate penalty
+[f,g,h] = st(r,s,k);
+
+
+function [f,g,h] = st(r,s,k)
+n = length(r);
+c = -n*(gammaln((k+1)/2) - gammaln(k/2) - .5*log(pi*s*k));
+f = c + .5*(k+1)*sum(log(1 + conj(r).*r/(s*k)));
+g = (k+1)*r./(s*k + conj(r).*r);
+h = (k+1)./(s*k + conj(r).*r);
|