bayes_dipole_trial

function [m,deviations] = bayes_dipole(a,v,sigma_v,sigma_m)

   a = forward coefficients (from forward model)
    v = data
    sigma_v = std deviation of data
   sigma_m = std deviation of moments (for prior)

    outputs
    m = dipole moments at maximum posterior likelihood
    deviations = std deviations
     deviations(1) = std deviation of potential errors 
          (should be ~= sigma_v)
     deviations(2) = std deviation of moments
          (should be ~= sigma_m)

0001 function [m,deviations] = bayes_dipole_trial(a,v,sigma_v,sigma_m)
0002 %function [m,deviations] = bayes_dipole(a,v,sigma_v,sigma_m)
0003 %
0004 %   a = forward coefficients (from forward model)
0005 %    v = data
0006 %    sigma_v = std deviation of data
0007 %   sigma_m = std deviation of moments (for prior)
0008 %
0009 %    outputs
0010 %    m = dipole moments at maximum posterior likelihood
0011 %    deviations = std deviations
0012 %     deviations(1) = std deviation of potential errors
0013 %          (should be ~= sigma_v)
0014 %     deviations(2) = std deviation of moments
0015 %          (should be ~= sigma_m)
0016 
0017 if ~(nargin == 4)
0018     error('improper parameter list')
0019 end;
0020 
0021 
0022 ls_mat = a'*a + ((sigma_v/sigma_m)^2  * eye(size(a,2)));
0023 
0024 inv_ls_mat = inv(ls_mat);
0025 lv_mat = a'*v;
0026 
0027 m = inv_ls_mat*lv_mat;
0028 
0029 error = a*m - v;
0030 
0031 deviations(1,:) = sqrt(sum(error.^2)/size(error,1));
0032 
0033 deviations(2,:) = sqrt(sum(m.^2)/size(m,1));