function DCM = J_dcm_model_estimate(DCM, doplot)
% This function conducts model fitting and estimation
% Input: DCM is the output of J_dcm_model_specify.m
%
% ---------- Frances Jin edited based on Karl Friston's script ------------
% October. 2021

rng('default')

if nargin < 2
    doplot = true;
end

% Unpack
M = DCM.M;

% multiple starting point model inversion using dcm fit
F = -Inf;
for i = 1:8
            
    % initialise parameters
    %------------------------------------------------------------------
    dP              = spm_sqrtm(M.pC)*randn(spm_length(M.pE))/4;
    DCM.M.pE        = spm_unvec(spm_vec(M.pE) + dP,M.pE);
    
    DCM.multistartcheck(i).pE = DCM.M.pE;
    % estimate candidate model
    %------------------------------------------------------------------
    dcm             = spm_dcm_fit(DCM, false);
    dcm             = dcm{1};
           
    % store model and update priors if evidence has increased
    %------------------------------------------------------------------
    if dcm.F > F
        F           = dcm.F;
        DCM         = dcm;
        M.pE        = DCM.Ep;
    end    
    DCM.multistartcheck(i).Ep = DCM.Ep;
    DCM.multistartcheck(i).F = F;
end

% generate predicted symptom scores using the posterior
%==========================================================================
Q = DCM.Ep;                                  % posterior expectation
U = DCM.xU;

% generate a synthetic timeseries
[sim_syx.y, sim_syx.s, sim_syx.x, sim_syx.U] = spm_nosology_gen(Q,M,U);
sim_syx.Y = DCM.xY;

% Store simulation results in DCM
DCM.y = sim_syx.y; % predicted data (symptoms)
DCM.s = sim_syx.s; % psychological states
DCM.x = sim_syx.x; % pathophysiological states
DCM.U = sim_syx.U; % estimated exogenous inputs
DCM.Y = sim_syx.Y; % empirical data

% Plot
%==========================================================================
if doplot
    J_dcm_plot(DCM);
end