%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Example MATLAB scripts for using BFDA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% --- Add pathes of the required MATLAB packages --- % BFDA, bspline, fdaM, mcmcdiag, PACE % Please REPLACE tool_dir by your BFDA package directory tool_dir = '../Software' ; addpath(genpath(cat(2, tool_dir, '/BFDA/Code'))) addpath(genpath(cat(2, tool_dir, '/bspline'))) addpath(genpath(cat(2, tool_dir, '/fdaM'))) addpath(genpath(cat(2, tool_dir, '/mcmcdiag'))) addpath(genpath(cat(2, tool_dir, '/PACErelease2.11'))) %% --- Set up parameters for simulation --- n = 30; % Number of functional data samples p = 40; % Number of pool-grid points, or evaluation-grid points s = sqrt(5); % Standard deviation of functional observations r = 2; % Signal-to-Noise ratio rho = 1/2; % Scale parameter in the Matern function nu = 3.5; % Order parameter in the Matern function pgrid = (0 : (pi/2)/(p-1) : (pi/2)); % Pool-grid dense = 0.6; % Proportion of the observation points of the pool-grid au = 0; bu = pi/2; % Observation domain of function data m = 20; % Number of working-grid points stat = 1; % Specify stationary data cgrid = 1; % Specify common observation grid %% Setup random seed stream = RandStream('twister','Seed', 2016); reset(stream); % set up a seed for simulation %% --- Analyzing stationary functional data with common grid --- % Generate simulation data from GP(3sin(4t), s^2 Matern_cor(d; rho, nu)) % with noises from N(0, (s/r)^2) GausFD_cgrid = sim_gfd(pgrid, n, s, r, nu, rho, dense, cgrid, stat); % Setup parameters for BFDA % Run with BHM param = setOptions_bfda('smethod', 'bhm', 'cgrid', 1, 'mat', 1, ... 'M', 10000, 'Burnin', 2000, 'w', 1, 'ws', 1); % Run with Bayesian Functional PCA % param = setOptions_bfda('smethod', 'bfpca', 'M', 50, 'Burnin', 20); % Run with standard Bayesian Gaussian Process model % param = setOptions_bfda('smethod', 'bgp', 'mat', 1, ... % 'M', 50, 'Burnin', 20); % Run with Cubic Smoothing Splines % param = setOptions_bfda('smethod', 'css', 'mat', 1, 'M', ... % 50, 'Burnin', 20, 'pace', 0); % call BFDA [out_cgrid, param] = ... BFDA(GausFD_cgrid.Xraw_cell, GausFD_cgrid.Tcell, param); %% --- Analyzing stationary functional data with uncommon grid --- GausFD_ucgrid = sim_gfd(pgrid, n, s, r, nu, rho, dense, 0, stat); param_uc = setOptions_bfda('smethod', 'bhm', 'cgrid', 0, 'mat', 1, 'M',... 10000, 'Burnin', 2000, 'pace', 1, 'ws', 0.1); [out_ucgrid, param_uc] = ... BFDA(GausFD_ucgrid.Xraw_cell, GausFD_ucgrid.Tcell, param_uc); %% --- Plot results with stationary data by BHM (Figure 1 and Figure 2) --- % Set up data for stationary-cgrid-ucgrid plot Xraw_cgrid = reshape(cell2mat(GausFD_cgrid.Xraw_cell), length(GausFD_cgrid.Xraw_cell{1}), n); Xtrue_cgrid = reshape(cell2mat(GausFD_cgrid.Xtrue_cell), length(GausFD_cgrid.Xtrue_cell{1}), n); T_cgrid = reshape(cell2mat(GausFD_cgrid.Tcell), length(GausFD_cgrid.Tcell{1}), n); Xfull_cgrid = Xraw_cgrid; mu_sample_cgrid = nanmean(Xfull_cgrid, 2); Xraw_ucgrid = reshape(cell2mat(GausFD_ucgrid.Xraw_cell), length(GausFD_ucgrid.Xraw_cell{1}), n); Xtrue_ucgrid = reshape(cell2mat(GausFD_ucgrid.Xtrue_cell), length(GausFD_ucgrid.Xtrue_cell{1}), n); T_ucgrid = reshape(cell2mat(GausFD_ucgrid.Tcell), length(GausFD_ucgrid.Tcell{1}), n); Xfull_ucgrid = NaN(p, n); % n by p data matrix with nan's for unobserved data for i = 1:n Idx = find(ismember(pgrid, GausFD_ucgrid.Tcell{i})); Xfull_ucgrid(Idx, i) = GausFD_ucgrid.Xraw_cell{i}; end mu_sample_ucgrid = nanmean(Xfull_ucgrid, 2); % Plot stationary functional data with common and uncommon grids (Figure 1) leglocation='northeast'; clims = [0, 1]; xlims = [0, pi/2]; ylims = [-10, 8]; h = figure('DefaultLegendFontSize',18,'DefaultLegendFontSizeMode','manual'); subplot_tight(2, 2, 1, [.1, .05]) idx = [1]; plot(T_cgrid(:, idx), Xraw_cgrid(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_cgrid(:, idx), Xtrue_cgrid(:, idx), 'r-.', pgrid, out_cgrid.Z(:, idx), 'b-', pgrid, out_cgrid.Z_CL(:, idx), 'b.', ... eval_grid, out_cgrid.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) [hleg1, hobj1] = legend('Sample', 'Truth', 'BHM', 'BHM 95% CI', 'Location', leglocation); idx = [6]; plot(T_cgrid(:, idx), Xraw_cgrid(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); plot(T_cgrid(:, idx), Xtrue_cgrid(:, idx), 'r-.', pgrid, out_cgrid.Z(:, idx), 'b-', pgrid, out_cgrid.Z_CL(:, idx), 'b.', ... eval_grid, out_cgrid.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off xlim(xlims) ylim(ylims) title('(a)'); %Functional data (common grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 2, [.1, .05]) idx = [1]; plot(T_ucgrid(:, idx), Xraw_ucgrid(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_ucgrid(:, idx), Xtrue_ucgrid(:, idx), 'r-.', pgrid, out_ucgrid.Z(:, idx), 'b-', pgrid, out_ucgrid.Z_CL(:, idx), 'b.', ... eval_grid, out_ucgrid.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) idx = [6]; plot(T_ucgrid(:, idx), Xraw_ucgrid(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); plot(T_ucgrid(:, idx), Xtrue_ucgrid(:, idx), 'r-.', pgrid, out_ucgrid.Z(:, idx), 'b-', pgrid, out_ucgrid.Z_CL(:, idx), 'b.', ... pgrid, out_ucgrid.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off xlim(xlims) ylim(ylims) title('(b)'); set(gca,'fontsize', 20) subplot_tight(2, 2, 3, [.1, .05]) plot(T_cgrid(:, 1), mu_sample_cgrid, 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_cgrid(:, 1), GausFD_cgrid.Mean_true, 'r-.', 'LineWidth', 2); hold on plot(T_cgrid(:, 1), out_cgrid.mu, 'b-', pgrid, out_cgrid.mu_CI, 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off title('(c)'); xlim(xlims) ylim([-5, 5]) set(gca,'fontsize', 20) subplot_tight(2, 2, 4, [.1, .05]) plot(pgrid, mu_sample_ucgrid, 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(pgrid, GausFD_ucgrid.Mean_true, 'r-.', 'LineWidth', 2); hold on plot(pgrid, out_ucgrid.mu, 'b-', pgrid, out_ucgrid.mu_CI, 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off title('(d)'); xlim(xlims) ylim([-5, 5]) set(gca,'fontsize', 20) set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/stat_sig')) % Plot covariance estimates (Figure 2) zlims = [-1, 8]; [xgrid,ygrid] = meshgrid(pgrid); % 3D plots of the covariance matrices h = figure(); subplot_tight(2, 2, 1, [.1,.05]) mesh(xgrid, ygrid, out_cgrid.Sigma) zlim(zlims) title('(a)'); % BHM covariance estimate(common grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 2, [.1,.05]) mesh(xgrid, ygrid, out_ucgrid.Sigma) zlim(zlims) title('(b)'); %BHM covariance estimate(uncommon grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 3, [.1,.05]) mesh(xgrid, ygrid, cov(Xfull_cgrid')) zlim(zlims) title('(c)'); %Sample covariance estimate(common grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 4, [.1,.05]) mesh(xgrid, ygrid, GausFD_ucgrid.Cov_true) zlim(zlims) title('(d)'); %True covariance(common/uncommon grids). set(gca,'fontsize', 20) set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/stat_cov')) %% --- Analyzing non-stationary functional data with common grid --- GausFD_cgrid_ns = sim_gfd(pgrid, n, s, r, nu, rho, dense, cgrid, 0); param_ns = setOptions_bfda('smethod', 'bhm', 'cgrid', 1, 'mat', 0, 'M',... 10000, 'Burnin', 2000, 'pace', 1, 'ws', 0.01); [out_cgrid_ns, param_ns] = ... BFDA(GausFD_cgrid_ns.Xraw_cell, GausFD_cgrid_ns.Tcell, param_ns); %% --- Analyzing non-stationary functional data with uncommon grid --- GausFD_ucgrid_ns = sim_gfd(pgrid, n, s, r, nu, rho, dense, 0, 0); param_uc_ns = setOptions_bfda('smethod', 'bhm', 'cgrid', 0, 'mat', 0, ... 'M', 10000, 'Burnin', 2000, 'pace', 1, 'ws', 0.01); [out_ucgrid_ns, param_uc_ns ] = ... BFDA(GausFD_ucgrid_ns.Xraw_cell, GausFD_ucgrid_ns.Tcell, param_uc_ns); %% --- Plot Results with non-stationary data by BHM (Figure 3 and Figure 4) --- % set up data for Nonstationary-cgrid-ucgrid plot Xraw_cgrid_ns = reshape(cell2mat(GausFD_cgrid_ns.Xraw_cell), length(GausFD_cgrid_ns.Xraw_cell{1}), n); Xtrue_cgrid_ns = reshape(cell2mat(GausFD_cgrid_ns.Xtrue_cell), length(GausFD_cgrid_ns.Xtrue_cell{1}), n); T_cgrid_ns = reshape(cell2mat(GausFD_cgrid_ns.Tcell), length(GausFD_cgrid_ns.Tcell{1}), n); Xfull_cgrid_ns = Xraw_cgrid_ns; mu_sample_cgrid_ns = nanmean(Xfull_cgrid_ns, 2); Xraw_ucgrid_ns = reshape(cell2mat(GausFD_ucgrid_ns.Xraw_cell), length(GausFD_ucgrid_ns.Xraw_cell{1}), n); Xtrue_ucgrid_ns = reshape(cell2mat(GausFD_ucgrid_ns.Xtrue_cell), length(GausFD_ucgrid_ns.Xtrue_cell{1}), n); T_ucgrid_ns = reshape(cell2mat(GausFD_ucgrid_ns.Tcell), length(GausFD_ucgrid_ns.Tcell{1}), n); Xfull_ucgrid_ns = NaN(p, n); % n by p data matrix with nan's for unobserved data for i = 1:n Idx = find(ismember(pgrid, GausFD_ucgrid_ns.Tcell{i})); Xfull_ucgrid_ns(Idx, i) = GausFD_ucgrid_ns.Xraw_cell{i}; end mu_sample_ucgrid_ns = nanmean(Xfull_ucgrid_ns, 2); % Plot results of nonstationary functional data by BHM (Figure 3) leglocation='southwest'; ylims = [-18, 18]; h = figure('DefaultLegendFontSize',18,'DefaultLegendFontSizeMode','manual'); subplot_tight(2, 2, 1, [.1, .05]) idx = [1]; plot(T_cgrid_ns(:, idx), Xraw_cgrid_ns(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_cgrid_ns(:, idx), Xtrue_cgrid_ns(:, idx), 'r-.', pgrid, out_cgrid_ns.Z(:, idx), 'b-', ... pgrid, out_cgrid_ns.Z_CL(:, idx), 'b.', ... pgrid, out_cgrid_ns.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) [hleg1, hobj1] = legend('Sample', 'Truth', 'BHM', 'BHM 95% CI', 'Location', leglocation); idx = [6]; plot(T_cgrid_ns(:, idx), Xraw_cgrid_ns(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); plot(T_cgrid_ns(:, idx), Xtrue_cgrid_ns(:, idx), 'r-.', pgrid, out_cgrid_ns.Z(:, idx), 'b-', ... pgrid, out_cgrid_ns.Z_CL(:, idx), 'b.', ... pgrid, out_cgrid_ns.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off xlim(xlims) ylim(ylims) title('(a)'); % Functional data (common grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 2, [.1, .05]) idx = [1]; plot(T_ucgrid_ns(:, idx), Xraw_ucgrid_ns(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_ucgrid_ns(:, idx), Xtrue_ucgrid_ns(:, idx), 'r-.', pgrid, out_ucgrid_ns.Z(:, idx), 'b-',... pgrid, out_ucgrid_ns.Z_CL(:, idx), 'b.', ... pgrid, out_ucgrid_ns.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) idx = [10]; plot(T_ucgrid_ns(:, idx), Xraw_ucgrid_ns(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); plot(T_ucgrid_ns(:, idx), Xtrue_ucgrid_ns(:, idx), 'r-.', pgrid, out_ucgrid_ns.Z(:, idx), 'b-', ... pgrid, out_ucgrid_ns.Z_CL(:, idx), 'b.', ... pgrid, out_ucgrid_ns.Z_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off xlim(xlims) ylim(ylims) title('(b)'); % Functional data (uncommon grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 3, [.1, .05]) plot(pgrid, mu_sample_cgrid_ns, 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(pgrid, GausFD_cgrid_ns.Mean_true, 'r-.', 'LineWidth', 2); hold on plot(pgrid, out_cgrid_ns.mu, 'b-', pgrid, out_cgrid_ns.mu_CI, 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off title('(c)'); % Means (common grids). xlim(xlims) ylim([-8, 8]) set(gca,'fontsize', 20) subplot_tight(2, 2, 4, [.1, .05]) plot(pgrid, mu_sample_ucgrid_ns, 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(pgrid, GausFD_ucgrid_ns.Mean_true, 'r-.', 'LineWidth', 2); hold on plot(pgrid, out_ucgrid_ns.mu, 'b-', pgrid, out_ucgrid_ns.mu_CI, 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off title('(d)'); % Means (uncommon grids). xlim(xlims) ylim([-8, 8]) set(gca,'fontsize', 20) set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/non_stat_sig')) % Plot covariance estimates (Figure 4) zlims = [-1, 35]; h = figure(); subplot_tight(2, 2, 1, [.1,.05]) mesh(xgrid, ygrid, out_cgrid_ns.Sigma) zlim(zlims) title('(a)'); % BHM covariance estimate (common grids). set(gca,'fontsize', 20) subplot_tight(2, 2, 2, [.1,.05]) mesh(xgrid, ygrid, out_ucgrid_ns.Sigma) zlim(zlims) title('(b)'); set(gca,'fontsize', 20) % BHM covariance estimate (uncommon grids). subplot_tight(2, 2, 3, [.1,.05]) mesh(xgrid, ygrid, cov(Xfull_cgrid_ns')) zlim(zlims) title('(c)'); set(gca,'fontsize', 20) % Sample covariance estimate (common grids). subplot_tight(2, 2, 4, [.1,.05]) mesh(xgrid, ygrid, GausFD_ucgrid_ns.Cov_true) zlim(zlims) title('(d)'); set(gca,'fontsize', 20) % True covariance (common/uncommon grids). set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/non_stat_cov')) %% --- Analyzing stationary functional data with random grids --- GausFD_rgrid = sim_gfd_rgrid(n, p, au, bu, s, r, nu, rho, stat); eval_grid = prctile(sort(unique(cell2mat(GausFD_rgrid.Tcell))), ... 1:(100/p):100); param_rgrid = setOptions_bfda('smethod', 'babf', 'cgrid', 0, 'mat', 1, ... 'M', 10000, 'Burnin', 2000, 'm', m, 'eval_grid', eval_grid, 'ws', 1); % call BFDA [out_rgrid, param_rgrid]= ... BFDA(GausFD_rgrid.Xraw_cell, GausFD_rgrid.Tcell, param_rgrid); %% --- Analyzing nonstationary functional data with random grids --- GausFD_rgrid_ns = sim_gfd_rgrid(n, p, au, bu, s, r, nu, rho, 0); eval_grid_ns = prctile(sort(unique(cell2mat(GausFD_rgrid_ns.Tcell))), ... 1:(100/p):100); param_rgrid_ns = setOptions_bfda('smethod', 'babf', 'cgrid', 0, 'mat', ... 0, 'M', 10000, 'Burnin', 2000, 'm', m, 'eval_grid', ... eval_grid_ns, 'ws', 0.05); % call BFDA [out_rgrid_ns, param_rgrid_ns] = ... BFDA(GausFD_rgrid_ns.Xraw_cell, GausFD_rgrid_ns.Tcell, param_rgrid_ns); %% --- Plot results with random grids by BABF (Figure 5 and Figure 6) --- % Setup data for random grids case Xraw_rgrid = reshape(cell2mat(GausFD_rgrid.Xraw_cell), length(GausFD_rgrid.Xraw_cell{1}), n); Xtrue_rgrid = reshape(cell2mat(GausFD_rgrid.Xtrue_cell), length(GausFD_rgrid.Xtrue_cell{1}), n); T_rgrid = reshape(cell2mat(GausFD_rgrid.Tcell), length(GausFD_rgrid.Tcell{1}), n); Xraw_rgrid_ns = reshape(cell2mat(GausFD_rgrid_ns.Xraw_cell), length(GausFD_rgrid_ns.Xraw_cell{1}), n); Xtrue_rgrid_ns = reshape(cell2mat(GausFD_rgrid_ns.Xtrue_cell), length(GausFD_rgrid_ns.Xtrue_cell{1}), n); T_rgrid_ns = reshape(cell2mat(GausFD_rgrid_ns.Tcell), length(GausFD_rgrid_ns.Tcell{1}), n); % plot stationary/nonstationary functional data with random grids (Figure 5) leglocation='southwest'; ylims = [-11, 11]; h = figure('DefaultLegendFontSize',18,'DefaultLegendFontSizeMode','manual'); subplot_tight(2, 2, 1, [.1,.05]) idx = [1]; plot(T_rgrid(:, idx), Xraw_rgrid(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_rgrid(:, idx), Xtrue_rgrid(:, idx), 'r-.', eval_grid, out_rgrid.Z_cgrid(:, idx), 'b-', ... eval_grid, out_rgrid.Z_cgrid_CL(:, idx), 'b.', ... eval_grid, out_rgrid.Z_cgrid_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) [hleg1, hobj1] = legend('Sample', 'Truth', 'BABF', 'BABF 95% CI', 'Location', leglocation); idx = [6]; plot(T_rgrid(:, idx), Xraw_rgrid(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); plot(T_rgrid(:, idx), Xtrue_rgrid(:, idx), 'r-.', eval_grid, out_rgrid.Z_cgrid(:, idx), 'b-', ... eval_grid, out_rgrid.Z_cgrid_CL(:, idx), 'b.', ... eval_grid, out_rgrid.Z_cgrid_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off xlim(xlims) ylim(ylims) title('(a)'); % Functional data (stationary). set(gca,'fontsize', 20) subplot_tight(2, 2, 2, [.1,.05]) idx = [1]; plot(T_rgrid_ns(:, idx), Xraw_rgrid_ns(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); hold on plot(T_rgrid_ns(:, idx), Xtrue_rgrid_ns(:, idx), 'r-.', eval_grid, out_rgrid_ns.Z_cgrid(:, idx), 'b-',... eval_grid, out_rgrid_ns.Z_cgrid_CL(:, idx), 'b.', ... eval_grid, out_rgrid_ns.Z_cgrid_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) idx = [10]; plot(T_rgrid_ns(:, idx), Xraw_rgrid_ns(:, idx), 'LineWidth', 2, 'Color', [0.75, 0.75, 0.75]); plot(T_rgrid_ns(:, idx), Xtrue_rgrid_ns(:, idx), 'r-.', eval_grid, out_rgrid_ns.Z_cgrid(:, idx), 'b-', ... eval_grid, out_rgrid_ns.Z_cgrid_CL(:, idx), 'b.', ... eval_grid, out_rgrid_ns.Z_cgrid_UL(:, idx), 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off xlim(xlims) ylim([-10, 13]) title('(b)'); % Functional data (nonstationary). set(gca,'fontsize', 20) subplot_tight(2, 2, 3, [.1,.05]) plot(GausFD_rgrid.Tcell{1}, GausFD_rgrid.Mean_true_cell{1}, 'r-.', 'LineWidth', 2); hold on plot(eval_grid, out_rgrid.mu_cgrid, 'b-', eval_grid, out_rgrid.mu_cgrid_CI, 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off title('(c)'); % Mean estimate (stationary). xlim(xlims) ylim([-8, 8]) set(gca,'fontsize', 20) subplot_tight(2, 2, 4, [.1,.05]) plot(GausFD_rgrid_ns.Tcell{1}, GausFD_rgrid_ns.Mean_true_cell{1}, 'r-.', 'LineWidth', 2); hold on plot(eval_grid, out_rgrid_ns.mu_cgrid, 'b-', eval_grid, out_rgrid_ns.mu_cgrid_CI, 'b.', 'LineWidth', 2, 'MarkerSize', 10) hold off title('(d)'); % Mean estimate (nonstationary). xlim(xlims) ylim([-8, 8]) set(gca,'fontsize', 20) set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/rand_sig')) % Plot covariance estimate (Figure 6) zlims = [-1, 25]; h = figure(); subplot_tight(2, 2, 1, [.1,.05]) mesh(xgrid, ygrid, out_rgrid.Sigma_cgrid) zlim([-1, 8]) title('(a)'); % BHM covariance estimate (stationary). set(gca,'fontsize', 20) subplot_tight(2, 2, 2, [.1,.05]) mesh(xgrid, ygrid, out_rgrid_ns.Sigma_cgrid) zlim(zlims) title('(b)'); % BHM covariance estimate (nonstationary). set(gca,'fontsize', 20) subplot_tight(2, 2, 3, [.1,.05]) mesh(xgrid, ygrid, GausFD_cgrid.Cov_true) zlim([-1, 8]) title('(c)'); % True covariance (stationary). set(gca,'fontsize', 20) subplot_tight(2, 2, 4, [.1,.05]) mesh(xgrid, ygrid, GausFD_cgrid_ns.Cov_true) zlim(zlims) title('(d)'); % True covariance (nonstationary). set(gca,'fontsize', 20) set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/Figures/rand_cov')) %% --- Calculate RMSE (root mean square error) --- display('RMSE of the estimated stationary covariance') rmse(out_cgrid.Sigma_SE, GausFD_cgrid.Cov_true) display('RMSE of the estimated functional data') Xtrue_mat = reshape(cell2mat(GausFD_cgrid.Xtrue_cell), ... size(GausFD_cgrid.Xtrue_cell{1}, 2), size(GausFD_cgrid.Xtrue_cell, 2)); rmse(out_cgrid.Z, Xtrue_mat) display('RMSE of the estimated non-stationary covariance') Ctrue_ns = cov_ns(pgrid, s, nu, rho); % calculate the true non-stationary covariance matrix rmse(out_cgrid_ns.Sigma_SE, Ctrue_ns) %% --- Save simulated data sets and BFDA results --- save('Simu_Data.mat', 'GausFD_cgrid', 'GausFD_ucgrid', ... 'GausFD_cgrid_ns', 'GausFD_ucgrid_ns', ... 'GausFD_rgrid', 'GausFD_rgrid_ns', 'p', ... 'eval_grid', 'eval_grid_ns', 'pgrid', 'n') save('Simu_Output.mat', 'out_cgrid', 'out_ucgrid', ... 'out_cgrid_ns', 'out_ucgrid_ns', ... 'out_rgrid', 'out_rgrid_ns') %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Example MATLAB scripts for functional linear regression using fdaM %%%% with functional data smoothed by BFDA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% --- Load functional data output from previous simulation analysis --- load('Simu_Data.mat'); load('Simu_Output.mat'); %% --- Set sample sizes, training data set, and test data set --- n = 30; % Number of functional curves p = 40; % Number of pooled grid points, or evaluated grid points au = 0; bu = pi/2; % Domain of t pgrid = (au : (bu)/(p-1) : bu); % Pooled grid trange = [au, bu]; sampind = sort(randsample(1:n,20,false)) ; samptest = find(~ismember(1:n, sampind)); n_train = length(sampind); n_test = length(samptest); cgrid = 0; Xtrue = zeros(p, n); Xraw = zeros(p, n); Xsmooth = zeros(p, n); if cgrid % Functional observations with common grids Xtrue = reshape(cell2mat(GausFD_cgrid.Xtrue_cell), p, n); Xsmooth = out_cgrid.Z(:, sampind); Xraw = reshape(cell2mat(GausFD_cgrid.Xraw_cell), p, n); else % Functional observations with random grids for i = 1:n xi = GausFD_rgrid.Xtrue_cell{i}; xrawi = GausFD_rgrid.Xraw_cell{i}; ti = GausFD_rgrid.Tcell{i}; %interpolate by cubic smoothing spline h = mean(diff(ti)); Xtrue(:, i) = csaps(ti, xi, 1/(1 + h^3/6), pgrid); Xraw(:, i) = csaps(ti, xrawi, 1/(1 + h^3/6), pgrid); zi = out_rgrid.Zt{i}; Xsmooth(:, i) = csaps(ti, zi, 1/(1 + h^3/6), pgrid); end % Xsmooth = out_rgrid.Z_cgrid; % Xsmooth(1, :) = Xsmooth(2, :) * 0.9; end Xtrain = Xsmooth(:, sampind); Xtest = Xsmooth(:, samptest); Xraw_train = Xraw(:, sampind); Xraw_test = Xraw(:, samptest); rmse(Xtrue, Xsmooth) rmse(Xtrue, Xraw) %% --- Generate response variables --- betamat = (pgrid') .^ 2 ; % --- Scalar respones deltat = pgrid(2)-pgrid(1); Avec_true = deltat.*(Xtrue'*betamat - ... 0.5.*(Xtrue(1, :)'.*betamat(1) + Xtrue(p,:)'.*betamat(p)) ); Avec = Avec_true + normrnd(0, 1, n, 1); Avec_train = Avec(sampind); Avec_test = Avec(samptest); Avec_train_true = Avec_true(sampind); Avec_test_true = Avec_true(samptest); % --- Functional responses ymat_true = Xtrue .* repmat(betamat, 1, n) ; ymat = ymat_true + normrnd(0, 1, p, n); ymat_train = ymat(:, sampind); ymat_test = ymat(:, samptest); ymat_train_true = ymat_true(:, sampind); ymat_test_true = ymat_true(:, samptest); %% --- Set up functional data structure xfd, yfd for fdaM --- xnbasis = 20; xbasis = create_bspline_basis(trange, xnbasis, 4); xfd_true = smooth_basis(pgrid, Xtrue, xbasis); xfd = smooth_basis(pgrid, Xtrain, xbasis); xfd_raw = smooth_basis(pgrid, Xraw_train, xbasis); [yfd_samp, df, gcv, beta, SSE, penmat, y2cMap, argvals, y] = ... smooth_basis(pgrid, ymat_train, xbasis); %% --- Set up the curvature penalty operator --- conbasis = create_constant_basis(trange); % create a constant basis wfd = fd([0, 1], conbasis); wfdcell = fd2cell(wfd); curvLfd = Lfd(2, wfdcell); % Set up xfdcell xfdcell = cell(1, 2); xfdcell{1} = fd(ones(1, n_train), conbasis); xfdcell{2} = xfd; xfd_raw_cell = cell(1, 2); xfd_raw_cell{1} = fd(ones(1, n_train), conbasis); xfd_raw_cell{2} = xfd_raw; % Set up betacell for scalar responses betafd0 = fd(0, conbasis); bnbasis = 10; betabasis = create_bspline_basis(trange, bnbasis, 4); betafd1 = fd(zeros(bnbasis, 1), betabasis); betacell_vecy = cell(1, 2); betacell_vecy{1} = fdPar(betafd0); betacell_vecy{2} = fdPar(betafd1, curvLfd, 0); % Set up betacell, yfd_par for functional responses betacell_fdy = cell(1, 2); betacell_fdy{1} = fdPar(betafd1, curvLfd, 0); betacell_fdy{2} = fdPar(betafd1, curvLfd, 0); yfd_par = fdPar(yfd_samp, curvLfd, 0); %% --- Compute cross-validated SSE's for a range of smoothing parameters --- %{ wt = ones(1, length(sampind)); lam = (0:0.1:1); nlam = length(lam); SSE_CV_vecy = zeros(nlam,1); SSE_CV_raw_vecy = zeros(nlam, 1); SSE_CV_fdy = zeros(nlam,1); SSE_CV_raw_fdy = zeros(nlam, 1); for ilam = 1:nlam; lambda_vecy = lam(ilam); betacelli_vecy = betacell_vecy; betacelli_vecy{2} = putlambda(betacell_vecy{2}, lambda_vecy); SSE_CV_vecy(ilam) = fRegress_CV(Avec_train, xfdcell, betacelli_vecy, wt); fprintf('Scalar responses, lambda %6.2f: SSE = %10.4f\n', ... lam(ilam), SSE_CV_vecy(ilam)); SSE_CV_raw_vecy(ilam) = fRegress_CV(Avec_train, xfd_raw_cell, ... betacelli_vecy, wt); fprintf('Scalar responses, lambda %6.2f: SSE = %10.4f\n', ... lam(ilam), SSE_CV_raw_vecy(ilam)); betacelli_fdy = betacell_fdy; betacelli_fdy{1} = putlambda(betacell_fdy{1}, lambda_vecy); betacelli_fdy{2} = putlambda(betacell_fdy{2}, lambda_vecy); yfd_par_i = putlambda(yfd_par, lambda_vecy); SSE_CV_fdy(ilam) = fRegress_CV(yfd_par_i, xfdcell, betacelli_fdy, wt); fprintf('Functional respones, lambda %6.2f: SSE = %10.4f\n', ... lam(ilam), SSE_CV_fdy(ilam)); SSE_CV_raw_fdy(ilam) = fRegress_CV(yfd_par_i, xfd_raw_cell, ... betacelli_fdy, wt); fprintf('Functional respones, lambda %6.2f: SSE = %10.4f\n', ... lam(ilam), SSE_CV_raw_fdy(ilam)); end %} %% --- Fit the linear model --- lambda = 0.1; wt = ones(1, length(sampind)); % --- Scalar responses betacell_vecy{2} = fdPar(betafd1, curvLfd, lambda); fRegressStruct_vecy = fRegress(Avec_train, xfdcell, betacell_vecy, wt); fRegressStruct_raw_vecy = ... fRegress(Avec_train, xfd_raw_cell, betacell_vecy, wt); % Get coefficients betaestcell_vecy = fRegressStruct_vecy.betahat; Avec_hat = fRegressStruct_vecy.yhat; intercept_vecy = getcoef(getfd(betaestcell_vecy{1})); disp(['Constant term = ',num2str(intercept_vecy)]) betaestcell_raw_vecy = fRegressStruct_raw_vecy.betahat; Avec_hat_raw = fRegressStruct_raw_vecy.yhat; intercept_raw_vecy = getcoef(getfd(betaestcell_raw_vecy{1})); disp(['Constant term = ',num2str(intercept_raw_vecy)]) display(['Scalar reponses:', 'fitted mse = ', ... num2str(mse(Avec_train_true, Avec_hat)), ... '; fitted mse_raw = ',num2str(mse(Avec_train_true, Avec_hat_raw))]) % Compute Rsquare covmat = cov([Avec_train, Avec_hat]); Rsqrd = covmat(1,2)^2/(covmat(1,1)*covmat(2,2)); disp(['R-squared = ',num2str(Rsqrd)]) covmat_raw = cov([Avec_train, Avec_hat_raw]); Rsqrd_raw = covmat_raw(1,2)^2/(covmat_raw(1,1)*covmat_raw(2,2)); disp(['raw R-squared = ',num2str(Rsqrd_raw)]) % Compute sigma resid_vecy = Avec_train - Avec_hat; SigmaE_vecy = mean(resid_vecy.^2); disp(['Scalar responses: SigmaE = ',num2str(SigmaE_vecy)]) resid_raw_vecy = Avec_train - Avec_hat_raw; SigmaE_raw_vecy = mean(resid_raw_vecy.^2); disp(['Scalar responses: Raw SigmaE = ',num2str(SigmaE_raw_vecy)]) % --- Functional responses betacell_fdy{1} = fdPar(betafd1, curvLfd, lambda); betacell_fdy{2} = fdPar(betafd1, curvLfd, lambda); yfd_par = fdPar(yfd_samp, curvLfd, lambda); fRegressStruct_fdy = fRegress(yfd_par, xfdcell, ... betacell_fdy, wt, y2cMap); fRegressStruct_raw_fdy = ... fRegress(yfd_par, xfd_raw_cell, betacell_fdy, wt, y2cMap); betaestcell_fdy = fRegressStruct_fdy.betahat; yfd_hat = fRegressStruct_fdy.yhat; intercept_fdy = eval_fd(pgrid, getfd(betaestcell_fdy{1})); betaestcell_raw_fdy = fRegressStruct_raw_fdy.betahat; yfd_hat_raw = fRegressStruct_raw_fdy.yhat; intercept_raw_fdy = eval_fd(pgrid, getfd(betaestcell_raw_fdy{1})); % MSE of fitted responses ymat_fitted = eval_fd(pgrid, yfd_hat); ymat_fitted_raw = eval_fd(pgrid, yfd_hat_raw); display(['mse = ', num2str(mse(ymat_train_true, ymat_fitted)), ... '; mse_raw = ',num2str(mse(ymat_train_true, ymat_fitted_raw))]) % Compute squared residual correlation resid_fdy = ymat_train_true - ymat_fitted; SigmaE_fdy = cov(resid_fdy'); resid_raw_fdy = ymat_train_true - ymat_fitted_raw; SigmaE_raw_fdy = cov(resid_raw_fdy'); %% --- Recompute the analysis to get confidence limits --- % --- Scalar responses stderrStruct_vecy = fRegress_stderr(fRegressStruct_vecy, ... eye(n_train), SigmaE_vecy); betastderrcell_vecy = stderrStruct_vecy.betastderr; stderrStruct_raw_vecy = ... fRegress_stderr(fRegressStruct_raw_vecy, eye(n_train), ... SigmaE_raw_vecy); betastderrcell_raw_vecy = stderrStruct_raw_vecy.betastderr; % Constant coefficient standard error: intercept_std_vecy = getcoef(betastderrcell_vecy{1}); intercept_ste_raw_vecy = getcoef(betastderrcell_raw_vecy{1}); % --- Functional responses stderrStruct_fdy = fRegress_stderr(fRegressStruct_fdy, y2cMap, ... SigmaE_fdy); % fixed a bug in fRegress_stderr.m at line 124: % bstderrfdj = data2fd(bstderrj, tfine, betabasisj); should be % bstderrfdj = data2fd(tfine, bstderrj, betabasisj); betastderrcell_fdy = stderrStruct_fdy.betastderr; stderrStruct_raw_fdy = ... fRegress_stderr(fRegressStruct_raw_fdy, y2cMap, SigmaE_raw_fdy); betastderrcell_raw_fdy = stderrStruct_raw_fdy.betastderr; % Coefficient standard error: intercept_std_fdy = eval_fd(pgrid, betastderrcell_fdy{1}); intercept_std_raw_fdy = eval_fd(pgrid, betastderrcell_raw_fdy{1}); %% --- Predict on test data --- %Set up xfd for test data xfd_test = smooth_basis(pgrid, Xtest, xbasis); xfd_raw_test = smooth_basis(pgrid, Xraw_test, xbasis); % --- Scalar responses xfdcell_test = cell(1, 2); xfdcell_test{1} = fd(ones(1, n_test), conbasis); xfdcell_test{2} = xfd_test; xfd_raw_test_cell = cell(1, 2); xfd_raw_test_cell{1} = fd(ones(1, n_test), conbasis); xfd_raw_test_cell{2} = xfd_raw_test; Avec_pred = fRegressPred(xfdcell_test, betaestcell_vecy); Avec_pred_raw = fRegressPred(xfd_raw_test_cell, betaestcell_raw_vecy); display(['Scalar responses predict mse = ', ... num2str(mse(Avec_test_true, Avec_pred)), ... '; Scalar responses with raw data predict mse_raw = ',... num2str(mse(Avec_test_true, Avec_pred_raw))]) % --- Functional responses ymat_test_pred = ... eval_fd(pgrid, fRegressPred(xfdcell_test, betaestcell_fdy, xbasis)); ymat_test_pred_raw = ... eval_fd(pgrid, fRegressPred(xfd_raw_test_cell, ... betaestcell_raw_fdy, xbasis)); display(['Functional response prediction mse = ', ... num2str(mse(ymat_test_true, ymat_test_pred)), ... '; Functional responses prediction with Raw data mse_raw = ',... num2str(mse(ymat_test_true, ymat_test_pred_raw))]) %% Generate Figures with functional linear regression analysis results (Figures 7-10) xlims = [0, pi/2]; % Plot one sample functional curve (Figure 7a) i=30; h = figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); subplot_tight(1, 1, 1, [.1, .1]) plot(GausFD_rgrid.Tcell{i}, GausFD_rgrid.Xraw_cell{i}, 'LineWidth', 3, 'Color', [0.75, 0.75, 0.75]) hold on plot(pgrid, Xsmooth(:, i), 'r-',... pgrid, Xraw(:, i), 'b--', pgrid, Xtrue(:, i), 'c:', ... 'LineWidth', 3, 'MarkerSize', 10) xlabel('t'); ylabel('x(t)'); title('(a)') %Example functional data xlim(xlims) [hleg, hobj] = legend('Raw data', 'BABF', 'CSS', 'Truth', 'Location', 'Best'); set(gca, 'fontsize', 26); hold off set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_sample_curve')) % plot the temperature coefficient function with scalar responses (Figure 7bc) h = figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); subplot_tight(2, 1, 1, [.1, .07]) plotbeta(betaestcell_vecy{2}, betastderrcell_vecy{2}) ylim([-2.5, 3]) title('(b)') %\beta(t) of Bayesian smoothed data hold on plot(pgrid, pgrid.^2, 'c:', 'Linewidth', 3) set(gca, 'fontsize', 26); hold off subplot_tight(2, 1, 2, [.1, .07]) plotbeta(betaestcell_raw_vecy{2}, betastderrcell_raw_vecy{2}) ylim([-2.5, 3]) title('(c)') % \beta(t) of cubic spline smoothed data hold on plot(pgrid, pgrid.^2, 'c:', 'Linewidth', 3) set(gca, 'fontsize', 26); hold off set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_vecy_beta')) % plot the fit with scalar repsonses (Figure 8a) h= figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); subplot_tight(1, 1, 1, [.1, .1]) plot(Avec_train_true, Avec_hat, 'rx', ... Avec_train_true, Avec_hat_raw, 'bo', ... Avec_train_true, Avec_train_true, 'c-', 'MarkerSize', 10, 'Linewidth', 3) ylabel(''); xlabel('Truth'); ylabel('Fitted'); title('(a)') % Fitted response values xlim([-5, 2]); [hleg, hobj] = legend('BABF', 'CSS', 'Location', 'NorthWest'); set(gca, 'fontsize', 26); set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_vecy_fitted')) % plot predict responses with scalar responses (8b) h=figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); %h=figure(); subplot_tight(1, 1, 1, [.1, .1]) plot(Avec_test_true, Avec_pred, 'rx', ... Avec_test_true, Avec_pred_raw, 'bo', ... Avec_test_true, Avec_test_true, 'c-', 'MarkerSize', 10, 'Linewidth', 3) ylabel(''); xlabel('Truth'); ylabel('Predicted'); title('(b)') %Predicted response values xlim([-5, 3]); [hleg, hobj] = legend('BABF', 'CSS', 'Location', 'NorthWest'); set(gca, 'fontsize', 26); set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_vecy_pred')) % plot intercept function with functional responses (Figure 9ac) h = figure('DefaultLegendFontSize',18,'DefaultLegendFontSizeMode','manual'); subplot_tight(2, 1, 1, [.1, .07]) plotbeta(betaestcell_fdy{1}, betastderrcell_fdy{1}) title('(a)') % Intercept \beta_0(t) of Bayesian smoothed data ylim([-1, 1]); set(gca, 'fontsize', 26); subplot_tight(2, 1, 2, [.1, .07]) plotbeta(betaestcell_raw_fdy{1}, betastderrcell_raw_fdy{1}) title('(c)') % Intercept \beta_0(t) of cubic spline smoothed data ylim([-1, 1]); set(gca, 'fontsize', 26); set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_fdy_intercept')) % plot the temperature coefficient function with functional responses (Figure 9 bd) h = figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); subplot_tight(2, 1, 1, [.1, .07]) plotbeta(betaestcell_fdy{2}, betastderrcell_fdy{2}) title('(b)') %\beta(t) of Bayesian smoothed data hold on plot(pgrid, pgrid.^2, 'c:', 'Linewidth',3) set(gca, 'fontsize', 26); ylim([-0.2, 3]); hold off subplot_tight(2, 1, 2, [.1, .07]) plotbeta(betaestcell_raw_fdy{2}, betastderrcell_raw_fdy{2}) title('(d)') % \beta(t) of cubic spline smoothed data hold on plot(pgrid, pgrid.^2, 'c:', 'Linewidth', 3) ylim([-0.2, 3]); set(gca, 'fontsize', 26); hold off set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_fdy_beta')) % Plot fitted output with functional responses (Figure 10a) h=figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); i=10; subplot_tight(1, 1, 1, [.1, .1]) plot(pgrid, ymat_fitted(:, i), 'r-', ... pgrid, ymat_fitted_raw(:, i), 'b--', ... pgrid, ymat_train_true(:, i), 'c:', 'LineWidth', 3) xlabel('t'); ylabel(''); xlim(xlims); title('(a)') % fitted functional responses [hleg, hobj] = legend('BABF', 'CSS', 'Truth', 'Location', 'SouthWest'); set(gca, 'fontsize', 26); set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_fdy_fitted')) % plot predict responses with functional responses (Figure 10b) h=figure('DefaultLegendFontSize',20,'DefaultLegendFontSizeMode','manual'); i=5; subplot_tight(1, 1, 1, [.1, .1]) plot(pgrid, ymat_test_pred(:, i), 'r-', ... pgrid, ymat_test_pred_raw(:, i), 'b--', ... pgrid, ymat_test_true(:, i), 'c:', 'LineWidth', 3) xlabel('t'); ylabel(''); xlim(xlims); title('(b)') % Predicted functional response set(gca, 'fontsize', 26); [hleg, hobj] = legend('BABF', 'CSS', 'Truth', 'Location', 'SouthWest'); set(h, 'PaperOrientation','landscape'); set(h, 'PaperUnits','normalized'); set(h, 'PaperPosition', [0 0 1.02 1]); print(h, '-dpdf', cat(2, pwd, '/reg_fdy_pred'))