Bayesian estimation of multifractal parameters: Input parameters
We detail the format and meaning of the optional outputs of the main functions of the toolbox.
Contents
Synthetic image data
The image used here for illustration is a realization of a random field with controlled multifractal properties: c1 is constant and equals c1=0.7, and c2 takes the values c2=-0.04 in side an elliptical region, and c2=-0.001 outside.
% add path to toolbox and load synthetic image clear all; clc; currentpath = cd;addpath([currentpath(1:end-5)]);addpath([currentpath(1:end-5) '/all_tools/']); load('MRWseries640_6.mat'); data=squeeze(data(:,:,5)); C1T=C1;C2T=C2; C1=[0.5 0.9];C2=[-0.05 0]; M2=MASK;M2(MASK==1)=C2T(1);M2(MASK==2)=C2T(2); % display data and theoretical values for c2 figure(1);clf;set(gcf,'position',[10 500 700 300]); subplot(121);imagesc(data);colormap bone; title('A synthetic image with two distinct zones of multifractality'); subplot(122); imagesc(M2); set(gca,'clim',C2); axis tight; title('Corresponding mask with prescribed values for c2');hold on; text('string',['c2=',num2str(C2T(1))],'units','normalized','pos',[0.5 0.8]); text('string',['c2=',num2str(C2T(2))],'color',[1 1 1]-1e-3,'units','normalized','pos',[0.5 0.45]);
Multifractal analysis parameters
4 mandatory inputs are required that determine how the multifractal analysis is performed: Choice of analysis wavelet; choice of global regularity parameter for wavelet leaders; lower and upper scale invariance cutoff. These parameters are inputs to "preproc_data_mvmfa.m".
Nwt = 2; % number of vanishing moments of wavelet (Daubechies) gamint = 0; % fractional integration parameter for wavelet leaders j1=1; % scaling range lower cutoff j2=2; % scaling range upper cutoff
Data and data handling parameters
For the data to be analyzed (in form of a single data vector / matrix / cube), the last dimensions is considered to be a "time" variable when applicable, all others are spatial. Whether the analysis is performed for a collection of time series or a collection of images (patches) is determined by the size (1 or 2 entries) of the input parameters sw and dsw described a bit further below. For instance, a data cube of size [n1 x n2 x n3] may be analyzed as a collection of n1 x n2 time series of length n3, or as a collection of n3 images of size n1 x n2.
The analysis window size is prescribed by the variabme "sw". It is either sw = [nt], where nt <= n3, for time-series, or sw = [nx ny], where nx<=n1 and ny<=n2, for images. The time series / images are analyzed for a covering of windows of size sw.
The overlap of the windows is determined by the mandatory data handling parameter "dsw", dsw=[dn1] for time series, or [dn1 dn2] for images. For instance, dsw = 0.5*[nx ny] where sw = [nx ny] yields analysis windows with 50% overlap in the x and in the y directions.
The data and the parameters "sw" and "dsw" are mandatory inputs to "preproc_data_mvmfa.m".
sw=[64 64]; % horizontal / vertical image patch size to be analyzed dsw=0.5*sw; % overlap of image patches (here, 50%) % Data cutting and computation of DWT and DFT [lxft,param_data] = preproc_data_mvmfa(data,sw,dsw,j1,j2,Nwt,gamint);
Optional parameters for the Bayesian model and estimators
The parameters of the Bayesian model and Gibbs sampler can be customized by passing on any of the following name-value pairs as additional inputs to the function "init_model_mvmfa.m", i.e. [param_model] = init_model_mvmfa(param_data, OPTIONAL_PARAMS); where "OPTIONAL_PARAMS" are the name-value pairs
OPTIONS:
- 'NIScale' : number of integral scales of the data (N+; default: 1). Decorrelation length of the data, given by LDecorr = sw / NIScale.
- 'Nbi' : number of burn-in samples of MCM chains (N+; default: 500)
- 'Nmc' : total number of samples of MCM chains (N+; default: 800)
- 'hpGMRF' : spatial / temporal / volume GMRF hyperparameters ([R+ R+ R+]; default: [8 16 4])
- 'MEMSAVE' : (0 or 1) - 1: low-memory computation of estimators without storing MCM chains (default: 1). Note: if you want to output the MCM chains for instection or post-processing, MEMSAVE must be set to 0.
- 'nanIP' : (0 or 1) - 1: output missing data parameters sampled from prior (default: 1). If the data contains missing samples (identified with NaNs), estimates for these data portions are nevertheless provided for these data portions by sampling from the prior distributions (i.e., no the data fit term).
ADVANCED and EXPERIMENTAL OPTIONS:
- 'rmMean' : use a mean-free model for c2 in the joint estimation for (c1,c2) (0 or 1; default: 0).
- 'aSAR' : estimate SAR hyperparameter (0 or 1; default: 1). NOTE: if you choose NOT to estimate the SAR hyperparameters, you can manually set it by overwriting the variable param_model.param.epi.
- 'PostVal' : output posterior value and MAP estimators estimate; ONLY for univariate / noninformative priors (0 or 1; default: 0). If PostVal=1, the second output MCM chain for c2) of est_bayes_mvmfa for univariate estimation (i.e., no regularization for c2 and c1) is a STRUCTURE with elements: Tc2 (posterior variance or MCM chains for c2, depending on MEMSAVE); c2MAP, c20MAP, c1MAP, c10MAP (MAP estimators for c2, c20, c1, c10); Pt (value of the log-posterior distribution corresponding to MAP estimates).
EXAMPLE of input name-vamue pairs: output also MAP estimators and posterior value, and MCM chains
[param_model] = init_model_mvmfa(param_data,'PostVal',1,'MEMSAVE',0); % Bayesian estimation for (c2,c1) using noninformative priors: MMSE and MAP estimators [estc2,TC2,estc20,TC20,estc1,TC1,estc10,TC10,estSAR] = est_bayes_mvmfa(lxft,param_model,[0 0]); % the MAP estimators and posterior value TC2 ec2{2}=estc2;ec1{2}=estc1; ec2{3}=TC2.c2MAP;ec1{3}=TC2.c1MAP; % linear fit (for comparison) [ec2{1},ec20{1},ec1{1},ec10{1}] = linfit_cp(param_data,0); % plot c2 L=sw(1);cc0=MASK(L/2:L:end,L/2:L:end);c1=zeros(size(cc0));c2=zeros(size(cc0)); for k=1:length(C1T);c1(cc0==k)=C1T(k);c2(cc0==k)=C2T(k);end; figure(2);clf; set(gcf,'position',[10 10 900 250]); subplot(141);imagesc(squeeze(c2));axis image; grid on; caxis(C2); title('c2 theo downsampled');xlabel('x');ylabel('y'); subplot(142);imagesc(squeeze(ec2{1}));axis image; grid on; caxis(C2); title('c2 linear regression');xlabel('x');ylabel('y'); subplot(143);imagesc(squeeze(ec2{2}));axis image; grid on; caxis(C2); title('c2 Bayesian IG MMSE');xlabel('x');ylabel('y'); subplot(144);imagesc(squeeze(ec2{3}));axis image; grid on; caxis(C2); title('c2 Bayesian IG MAP');xlabel('x');ylabel('y'); axes('Units','Normal');h = title('c2 - Patch(x,y)');set(gca,'visible','off');set(h,'visible','on');
TC2 =
Tc2: [361x800 double]
c2MAP: [19x19 double]
c20MAP: [19x19 double]
c1MAP: [19x19 double]
c10MAP: [19x19 double]
Pt: [19x19 double]
the log-posterior value an be used as a model selection criterion. For example, the synthetic image used in the above example has one integral scale (within homogeneous zones). A model with a much larger number of integral scales should yield a lower log-posterior value:
[param_model] = init_model_mvmfa(param_data,'NIScale',1,'PostVal',1); [estc2_1,ST1] = est_bayes_mvmfa(lxft,param_model); [param_model] = init_model_mvmfa(param_data,'NIScale',10,'PostVal',1); [estc2_2,ST2] = est_bayes_mvmfa(lxft,param_model); diff_post1_post2=sum(ST1.Pt(:)-ST2.Pt(:))
diff_post1_post2 = 245.1840