xiscalina / optimizer Goto Github PK

Optimizer tools written in Matlab. Including MCMC, DEMC, Ensemble Kalman filter, Approximate Bayesian Computing-Population Monte Carlo, and modeling averaging methods. Most of the algorithms are based on class CEE 290 taught by Prof. Jasper A. Vrugt @ UCI.

All input and output of functions are single structures with multiple fields. Same field names were used across different functions and can be piped together for testing purpose.

Differential Evolution / Differential Evolution Markov Chain. In the case of 'DE', this function also plot a figure for the objective function of all generations.

• func : evaluation (objective) function.
• bound : 2 X n matrix, where n is the parameter dimension. Lower (first row) and upper(second row) bounds of parameter.
• size : size of population.
• life : maximum generation of population.
• funcPrior :(optional) prior distribution function to generate first generation. Default is uninformative prior within bounds.
• type :(optional) 'DE' or 'DEMC'. The former seek the global minimum, the later maps the posterior. Default is 'DEMC'.
• greedy :(optional) logical scalar, faster convergence but less stable. Breaks detailed balance and should not be used with 'DEMC'.
• logFlag :(optional) logical scalar, use log liklihood or not. Default is 'false'.
• mutation :(optional) scalar, mutation rate.
• CR :(optional) scalar, crossover rate.
• simplex :(optional) 1 X n logical maxtrix. 'true' specifies the parameters that should be positive and sum up to 1.
• name :(optional) name of output figure. Default 'DEMCtest'.

• chain : population at each generation.
• obj : evaluation (objective) function value at each generation.
• best : best individual at last generation.

Markov Chain Monte Carlo with Random-Walk Metropolis.

One of 'loc', 'funcPrior', 'bound'(in descending priority) must be specified.

• life : length of Markov Chain.
• func : liklihood function. Calculates the liklihood from parameter.
• funcProp : (optional) proposal distribution function. Make new proposal
• according to parameter. Default is normal distribution with sigma = 1.
• loc : (optional) initial parameter.
• funcPrior : (optional) prior function to generate initial parameter.
• bound : (optional) 2 X n matrix, where n is the parameter dimension.
• Lower(first row) and upper(second row) bounds of parameter.

• chain : m X n matrix, where m = length and n = parameter dimension.
• time : m X 1 matrix, time stamp of each step.

Ensemble Kalman filter. Updated forecast and analysis states using observation and ensemble Kalman filter.

• R : measurement error (sigma).
• Fstates : forecast states. 1 X n matrix, n = number of state
• Astates : analysis states. 1 X n matrix.
• obs : observation states. Scaler.
• func : model operator, from current state to next state.

The same structue with updated Fstates and Astates, as well as estimation of model error (C).

Multi-model averaging function. Supported methods:

• EWA : Equal Weights Averaging
• BGA : Bates-Granger Averaging
• AICA : Akaike's Information Criterion Averaging
• BICA : Bayes Information Criterion Averaging
• GRA : Granger-Ramanathan Averaging
• BMA : Bayesian Model Averaging
• MMA : Mallows Model Averaging
• MMAd : Mallows Model Averaging (lie on simplex)

• Xcali : m1 X n model prediction matrix for calibration period. m1 = time
• length of calibration period, n = model number.
• Ycali : m1 X 1 observation matrix for calibration period.
• Xeval : m2 X n model prediction matrix for evaluation period.
• Yeval : m2 X 1 observation matrix for evaluation period.
• p : 1 X n matrix of model parameter numbers.
• method : abbreviation of method.

• weight : model weights.
• RMSEcali : RMSE of calibration period.
• RMSEeval : RMSE of evaluation period.
• chain : (BMA, MMA and MMAd) model weights in all iterations.
• obj : (BMA, MMA and MMAd) objective function value for all iterations.
• sigma : (BMA) optimized model sigma.

Approximate Bayesian Computing-Population Monte Carlo method.

• obs : 1 X d, d = dimension of target
• bound : 2 X n matrix, where n is the parameter dimension. Lower
• (first row) and upper(second row) bounds of parameter.
• funcPrior : prior distribution function to generate parameter.
• funcSummar : function to model target from parameter.
• funcDist : function to calculate distance between modeled and observed
• target.
• epsl : 1 X i epsilon matrix, where i is the number of generations.
• size : population size

The same strucuture with extra fields.

• inds : individuals (parameters) for all generations.
• summar : modeled target at all generations.
• weight : weight matrix for all generations.
• C : Covariance matrix for all geneations.
• acRate : overall acceptance rate.

For DEMC and MCMC:

%% Finding global minimum with DE

% Test function: 2d Ackley's function
pop.func = @(x) -20*exp(-.2*sqrt(.5*sum(x.^2)))...
    -exp(.5*sum(cos(2*pi*x)))+exp(1)+20;
pop.bound = [-5,-5;5,5];

% Setup the DE settings
pop.size = 50; pop.life = 200; pop.type = 'DE';
result = DEMC(pop); 
disp(['Global minimum from DE: ' num2str(result.best)])


%% Mapping posterior with MCMC and DEMC

% Test function: wide spread mixture of two normal distributions: 1/3 of
% N(-5,1) and 2/3 of N(5,1). We expect a mean of 1.67 and a std of 4.82
pop.func = @(x) 1/3*normpdf(x,-5,1)+2/3*normpdf(x,5,1);
pop.bound = [-100;100]+1.67;

% Using DEMC. It should be able to sample both peaks.
pop.type = 'DEMC'; pop.life = 1500;
result = DEMC(pop); burnt = result.chain(pop.life/2:end,:);
disp(['Mean estimated by DEMC: ' num2str(mean(burnt(:)))])
disp(['Std estimated by DEMC: ' num2str(std(burnt(:)))])

% Using MCMC. It should only be able to sample one peak (-5 or 5).
pop.life = 10000; pop.loc = 0;
result = MCMC(pop); burnt = result.chain(pop.life/2:end,:);
disp(['Mean estimated by MCMC: ' num2str(mean(burnt))])
disp(['Std estimated by MCMC: ' num2str(std(burnt))])