function [L,gam,momx]=loglios(th,params,Hk,k,scl) % [L,gam,momx]=LOGLIOS(th,params,Hk,k,scl) % % Calculates the full negative logarithmic likelihood and its % derivatives, i.e. minus LKOS and minus GAMMAKOS averaged over % wavenumber space. This is the function that we need to MINIMIZE! % % INPUT: % % th The five-parameter vector argument [scaled]: % th(1)=D Isotropic flexural rigidity % th(2)=f2 The sub-surface to surface initial loading ratio % th(3)=s2 The first Matern parameter, aka sigma^2 % th(4)=nu The second Matern parameter % th(5)=rho The third Matern parameter % params A structure with AT LEAST these constants that are known: % DEL surface and subsurface density contrast [kg/m^3] % g gravitational acceleration [m/s^2] % blurs 0 Don't blur likelihood using the Fejer window % N Blur likelihood using the Fejer window [default: N=2] % kiso wavenumber beyond which we are not considering the likelihood % Hk A [prod(params.NyNx)*2]-column vector of complex Fourier-domain observations % k the wavenumbers at which these are being evaluated [1/m] % scl The vector with any scalings applied to the parameter vector % momx Moments of the quadratic piece Xk over all relevant wavenumbers % % OUTPUT: % % L The loglihood, averaged over all wavenumbers % gam The score, averaged over all wavenumbers % % SEE ALSO: % % FISHERKOS, which should be incorporated at a later stage % % Last modified by fjsimons-at-alum.mit.edu, 10/22/2014 % Default scaling is none defval('scl',ones(size(th))) % Scale up the parameter vector for the proper likelihood and score th=th.*scl; % Here I build the protection that the flexural rigidity, % subsurface-to-surface ration, and the three Matern parameters should be % positive. I mirror them up! Thereby messing with the iteration path, % but hey. It means we can use FMINUNC also. th([1 2 3 4 5])=abs(th([1 2 3 4 5])); % Filter, perhaps [Lk,Xk]=Lkos(k,th,params,Hk); if any(~isnan(params.kiso)) Lk(k>params.kiso)=NaN; Xk(k>params.kiso)=NaN; end % Note: should we restrict this to the upper halfplane? or will mean do % Get the likelihood at the individual wavenumbers; average L=-nanmean(Lk); if isnan(L) % Attempt to reset L=1e100; end if nargout==3 % Extract the moments we'll be needing for evaluation later df=2; % First should be close to df/2, second close to df/2, third is like % the second except formally from a distribution that should be normal % with mean df/2 and variance blabla/K; the last parameter is "magic" momx=[nanmean(Xk) nanvar(Xk) nanmean([Xk-df/2].^2)]; end % I say, time to extract heskosl and hes here also? % Get the scores at the individual wavenumbers; average switch params.blurs case {0,1} gam=-nanmean(gammakos(k,th,params,Hk)); % The correct gradient is too heterogeneous to be good so scale gam=gam.*scl; otherwise gam=NaN; end % Print the trajectory, seems like one element at a time gets changed % disp(sprintf('Current theta: %8.3g %8.3g %8.3g %8.3g %8.3g',th))