diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc index 457e57b8b..f860cbe4f 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.cc +++ b/mex/sources/estimation/LogLikelihoodSubSample.cc @@ -176,9 +176,6 @@ LogLikelihoodSubSample::updateParams(const Vector &estParams, Vector &deepParams info = 1; } // end switch } // end found -#ifdef DEBUG - mexPrintf("End of Setting of HQ params\n"); -#endif } //end for }; diff --git a/mex/sources/estimation/RandSampler.hh b/mex/sources/estimation/RandSampler.hh index 16a08e427..5cd245d99 100644 --- a/mex/sources/estimation/RandSampler.hh +++ b/mex/sources/estimation/RandSampler.hh @@ -31,22 +31,18 @@ class RandSampler { public: - RandSampler(){ - }; virtual ~RandSampler(){}; - virtual double compute(Vector &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState, - Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, - size_t presampleStart, int &info, const size_t nMHruns = 0, const Matrix &Jscale, - const Matrix &D, const LogPosteriorDensity &logPosteriorDensity, const Prior &drawDistribution) = 0; - + virtual double compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState, + Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, + const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Jscale, + LogPosteriorDensity &logPosteriorDensity, Prior &drawDistribution, + EstimatedParametersDescription &epd) = 0; /** * draw = Mean + randn(1,n) * Sigma_upper_chol; * */ virtual void randMultiVar(Prior &distribution, Vector &draw, const Vector &mean, const Matrix &Scale, const size_t n); - virtual void saveDraws() = 0; - }; #endif // !defined(5D7E5E52_2A4F_4f98_9A5C_A7FD8C278E0A__INCLUDED_) diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.cc b/mex/sources/estimation/RandomWalkMetropolisHastings.cc index b92bb268b..5721bcf42 100644 --- a/mex/sources/estimation/RandomWalkMetropolisHastings.cc +++ b/mex/sources/estimation/RandomWalkMetropolisHastings.cc @@ -20,29 +20,43 @@ #include "RandomWalkMetropolisHastings.hh" double -RandomWalkMetropolisHastings::compute(Vector &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState, - Vector &estParams2, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, - const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Jscale, - const Matrix &D, LogPosteriorDensity &lpd, Prior &drawDistribution) +RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState, + Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, + const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Dscale, + LogPosteriorDensity &lpd, Prior &drawDistribution, EstimatedParametersDescription &epd) { - double logpost, newLogpost; - size_t accepted = 0; - parDraw = estParams2; - blas::gemm("N", "N", 1.0, D, Jscale, 1.0, Dscale); + bool overbound; + double newLogpost, logpost; + size_t count, accepted = 0; + parDraw = estParams; + logpost = - lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info); for (size_t run = 0; run < nMHruns; ++run) { + overbound=false; randMultiVar(drawDistribution, newParDraw, parDraw, Dscale, parDraw.getSize()); - try + for (count=0;count epd.estParams[count].upper_bound ); + if (overbound) + { + newLogpost = -INFINITY; + break; + } } - catch(...) + if (!overbound) { - newLogpost = -INFINITY; + try + { + newLogpost = - lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info); + } + catch(...) + { + newLogpost = -INFINITY; + } } if ((newLogpost > -INFINITY) && log(uniform.drand()) < newLogpost-logpost) { - mat::get_col(mhParams, run) = newParDraw; + mat::get_row(mhParams, run) = newParDraw; parDraw = newParDraw; mhLogPostDens(run) = newLogpost; logpost = newLogpost; @@ -50,16 +64,10 @@ RandomWalkMetropolisHastings::compute(Vector &mhLogPostDens, MatrixView &mhParam } else { - mat::get_col(mhParams, run) = parDraw; + mat::get_row(mhParams, run) = parDraw; mhLogPostDens(run) = logpost; } } - return accepted/nMHruns; -} - -void -RandomWalkMetropolisHastings::saveDraws(const std::string &modName, const std::string &suffix, const MatrixView &Draws, const size_t block) -{ - + return (double) accepted/nMHruns; } diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.hh b/mex/sources/estimation/RandomWalkMetropolisHastings.hh index 7507f4ef8..bc64c98e8 100644 --- a/mex/sources/estimation/RandomWalkMetropolisHastings.hh +++ b/mex/sources/estimation/RandomWalkMetropolisHastings.hh @@ -28,27 +28,25 @@ #include "RandSampler.hh" -class RandomWalkMetropolisHastings : public RandSampler { +class RandomWalkMetropolisHastings : public RandSampler +{ private: UniformPrior uniform; - Matrix Dscale; Vector parDraw, newParDraw; public: RandomWalkMetropolisHastings(size_t size) : uniform(0.0, 0.0, 0.0, 1.0, 0.0, 1.0), - Dscale(size), parDraw(size), newParDraw(size) + parDraw(size), newParDraw(size) { }; virtual ~RandomWalkMetropolisHastings(){}; - double compute(Vector &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState, - Vector &estParams2, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, - const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Jscale, - const Matrix &D, LogPosteriorDensity &logPosteriorDensity, Prior &drawDistribution); -// void randMultiVar(const Prior &distribution, Vector &draw, const Matrix &Scale, const Vector &mean, const size_t n = 0); - void saveDraws(const std::string &modName, const std::string &suffix, const MatrixView &Draws, const size_t block); - + virtual double compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState, + Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, + const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Jscale, + LogPosteriorDensity &logPosteriorDensity, Prior &drawDistribution, + EstimatedParametersDescription &epd); }; #endif // !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_) diff --git a/mex/sources/estimation/libmat/Vector.hh b/mex/sources/estimation/libmat/Vector.hh index cf42438ce..cfedb52e7 100644 --- a/mex/sources/estimation/libmat/Vector.hh +++ b/mex/sources/estimation/libmat/Vector.hh @@ -227,6 +227,27 @@ namespace vec } return nrm; } + + // Computes the sum, min and max of a vector + // and returns double mean=sum/size + template + double + meanSumMinMax(double &sum, double &min, double &max, const Vec &v) + { + sum = 0; + min=max=v(0); + const double *p = v.getData(); + while (p < v.getData() + v.getSize() * v.getStride()) + { + if ((*p) > max) + max = (*p); + if ((*p) < min) + min = (*p); + sum+=*p; + p += v.getStride(); + } + return sum/v.getSize(); + } } // End of namespace #endif diff --git a/mex/sources/estimation/logMHMCMCposterior.cc b/mex/sources/estimation/logMHMCMCposterior.cc new file mode 100644 index 000000000..94d626e47 --- /dev/null +++ b/mex/sources/estimation/logMHMCMCposterior.cc @@ -0,0 +1,451 @@ +/* + * Copyright (C) 2010 Dynare Team + * + * This file is part of Dynare. + * + * Dynare is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Dynare is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Dynare. If not, see . + */ + +#include +#include +#include +#include + +#include "Vector.hh" +#include "Matrix.hh" +#include "LogPosteriorDensity.hh" +#include "RandomWalkMetropolisHastings.hh" + +#include "mex.h" +#include "mat.h" + +#if defined(_WIN32) || defined(__CYGWIN32__) || defined(WINDOWS) +# define DIRECTORY_SEPARATOR "\\" +#else +# define DIRECTORY_SEPARATOR "/" +#endif + +void +fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType type, + std::vector &estParamsInfo) +{ + // execute once only + static const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_"); + static const mxArray *bayestopt_ubp = mxGetField(bayestopt_, 0, "ub"); // upper bound + static const mxArray *bayestopt_lbp = mxGetField(bayestopt_, 0, "lb"); // lower bound + static const mxArray *bayestopt_p1p = mxGetField(bayestopt_, 0, "p1"); // prior mean + static const mxArray *bayestopt_p2p = mxGetField(bayestopt_, 0, "p2"); // prior standard deviation + static const mxArray *bayestopt_p3p = mxGetField(bayestopt_, 0, "p3"); // lower bound + static const mxArray *bayestopt_p4p = mxGetField(bayestopt_, 0, "p4"); // upper bound + static const mxArray *bayestopt_p6p = mxGetField(bayestopt_, 0, "p6"); // first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution). + static const mxArray *bayestopt_p7p = mxGetField(bayestopt_, 0, "p7"); // second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution). + static const mxArray *bayestopt_jscalep = mxGetField(bayestopt_, 0, "jscale"); // MCMC jump scale + + static const size_t bayestopt_size = mxGetM(bayestopt_); + static const VectorConstView bayestopt_ub(mxGetPr(bayestopt_ubp), bayestopt_size, 1); + static const VectorConstView bayestopt_lb(mxGetPr(bayestopt_lbp), bayestopt_size, 1); + static const VectorConstView bayestopt_p1(mxGetPr(bayestopt_p1p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p1"); + static const VectorConstView bayestopt_p2(mxGetPr(bayestopt_p2p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p2"); + static const VectorConstView bayestopt_p3(mxGetPr(bayestopt_p3p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p3"); + static const VectorConstView bayestopt_p4(mxGetPr(bayestopt_p4p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p4"); + static const VectorConstView bayestopt_p6(mxGetPr(bayestopt_p6p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p6"); + static const VectorConstView bayestopt_p7(mxGetPr(bayestopt_p7p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p7"); + static const VectorConstView bayestopt_jscale(mxGetPr(bayestopt_jscalep), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "jscale"); + + // loop processsing + size_t m = mxGetM(estim_params_info), n = mxGetN(estim_params_info); + MatrixConstView epi(mxGetPr(estim_params_info), m, n, m); + size_t bayestopt_count = estParamsInfo.size(); + + for (size_t i = 0; i < m; i++) + { + size_t col = 0; + size_t id1 = (size_t) epi(i, col++) - 1; + size_t id2 = 0; + if (type == EstimatedParameter::shock_Corr + || type == EstimatedParameter::measureErr_Corr) + id2 = (size_t) epi(i, col++) - 1; + col++; // Skip init_val #2 or #3 + double par_low_bound = bayestopt_lb(bayestopt_count); col++; //#3 epi(i, col++); + double par_up_bound = bayestopt_ub(bayestopt_count); col++; //#4 epi(i, col++); + Prior::pShape shape = (Prior::pShape) epi(i, col++); + double mean = epi(i, col++); + double std = epi(i, col++); + double low_bound = bayestopt_p3(bayestopt_count); + double up_bound = bayestopt_p4(bayestopt_count); + double fhp = bayestopt_p6(bayestopt_count); // double p3 = epi(i, col++); + double shp = bayestopt_p7(bayestopt_count); // double p4 = epi(i, col++); + + Prior *p = Prior::constructPrior(shape, mean, std, low_bound, up_bound, fhp, shp); //1.0,INFINITY);//p3, p4); + + // Only one subsample + std::vector subSampleIDs; + subSampleIDs.push_back(0); + estParamsInfo.push_back(EstimatedParameter(type, id1, id2, subSampleIDs, + par_low_bound, par_up_bound, p)); + bayestopt_count++; + } +} + +size_t +sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh, + Matrix &steadyState, Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, + const size_t presampleStart, int &info, const VectorConstView &nruns, const size_t fblock, const size_t nBlocks, + const Matrix &Jscale, Prior &drawDistribution, EstimatedParametersDescription &epd, + const std::string &resultsFileStem, const size_t console_mode) +{ + enum{iMin,iMax}; + std::vector OpenOldFile(nBlocks, 1); + size_t jloop = 0, irun, j; // counters + double dsum, dmax, dmin, sux = 0, jsux=0; + std::string mhFName; + std::stringstream ssFName; + MATFile *drawmat; // MCMC draws output file pointer + FILE *fidlog; // log file + int matfStatus; + size_t npar = estParams.getSize(); + Matrix MinMax(npar,2); + const mxArray *myinputs = mexGetVariablePtr("caller", "myinputs"); + const mxArray *InitSizeArrayPtr = mxGetField(myinputs, 0, "InitSizeArray"); + const VectorConstView InitSizeArrayVw(mxGetPr(InitSizeArrayPtr), nBlocks, 1); + Vector InitSizeArray(InitSizeArrayVw.getSize()); + InitSizeArray = InitSizeArrayVw; + const mxArray *flinePtr = mxGetField(myinputs, 0, "fline"); + const VectorConstView fline(mxGetPr(flinePtr), nBlocks, 1); + + const mxArray *NewFileArrayPtr = mxGetField(myinputs, 0, "NewFile"); + VectorView NewFileVw(mxGetPr(NewFileArrayPtr), nBlocks, 1); + Vector NewFile(NewFileVw.getSize()); + NewFile = NewFileVw; + + const mxArray *MAX_nrunsPtr = mxGetField(myinputs, 0, "MAX_nruns"); + const size_t MAX_nruns = (size_t) mxGetScalar(MAX_nrunsPtr); + + const mxArray *record = mxGetField(myinputs, 0, "record"); + const mxArray *AcceptationRatesPtr = mxGetField(record, 0, "AcceptationRates"); + VectorView AcceptationRates(mxGetPr(AcceptationRatesPtr), nBlocks, 1); + + mxArray *mxLastParametersPtr= mxGetField(record, 0, "LastParameters"); + MatrixView LastParameters(mxGetPr(mxLastParametersPtr), nBlocks, npar, nBlocks); + + mxArray *mxLastLogLikPtr = mxGetField(record, 0, "LastLogLiK"); + VectorView LastLogLiK(mxGetPr(mxLastLogLikPtr), nBlocks, 1); + + mxArray *mxMhLogPostDensPtr = 0; + mxArray *mxMhParamDrawsPtr = 0; + size_t currInitSizeArray = 0; + + mexPrintf("\n Starting MH Block Loop\n\n"); + + for (size_t b = fblock; b <= nBlocks; ++b) + { + jloop = jloop+1; + + sux = 0.0; + jsux = 0; + irun = (size_t) fline(b-1); + j = 0;//1; + while (j < nruns(b-1)) + { + if (currInitSizeArray != (size_t) InitSizeArray(b-1)) + { + // new or different size result arrays/matrices + currInitSizeArray = (size_t) InitSizeArray(b-1); + if (mxMhLogPostDensPtr) + mxDestroyArray(mxMhLogPostDensPtr); // log post density array + mxMhLogPostDensPtr = mxCreateDoubleMatrix(currInitSizeArray, 1, mxREAL); + if (mxMhParamDrawsPtr) + mxDestroyArray(mxMhParamDrawsPtr); // accepted MCMC MH draws + mxMhParamDrawsPtr = mxCreateDoubleMatrix( currInitSizeArray, npar, mxREAL); + } + VectorView mhLogPostDens(mxGetPr(mxMhLogPostDensPtr), currInitSizeArray, (size_t) 1); + MatrixView mhParamDraws(mxGetPr(mxMhParamDrawsPtr), currInitSizeArray, npar, currInitSizeArray); + jsux= rwmh.compute(mhLogPostDens, mhParamDraws, steadyState, estParams, deepParams, data, Q, H, + presampleStart, info, currInitSizeArray, Jscale, lpd, drawDistribution, epd); + sux+=jsux*currInitSizeArray; + j += currInitSizeArray; //j=j+1; + irun += currInitSizeArray; + + if(console_mode) + mexPrintf(" MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acceptance rate: %3.f \b%%\r", b, nBlocks, 100 * j/nruns(b-1), 100 * sux / j); + // % Now I save the simulations + // save draw 2 mat file ([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); + ssFName.clear(); + ssFName.str(""); + ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << resultsFileStem << "_mh" << (size_t) NewFile(b-1) << "_blck" << b << ".mat" ; + mhFName = ssFName.str(); + drawmat = matOpen(mhFName.c_str(), "w"); + if (drawmat==0) + { + mexPrintf("Error in MH: Can not open draws Mat file for writing: %s \n", mhFName.c_str()); + exit(1); + } + matfStatus = matPutVariable(drawmat, "x2", mxMhParamDrawsPtr); + if (matfStatus) + { + mexPrintf("Error in MH: Can not use draws Mat file for writing: %s \n", mhFName.c_str()); + exit(1); + } + matfStatus = matPutVariable(drawmat, "logpo2", mxMhLogPostDensPtr); + if (matfStatus) + { + mexPrintf("Error in MH: Can not usee draws Mat file for writing: %s \n", mhFName.c_str()); + exit(1); + } + matClose(drawmat); + + // save log to fidlog = fopen([MhDirectoryName '/metropolis.log'],'a'); + ssFName.str(""); + ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << "metropolis.log" ; + mhFName = ssFName.str(); + fidlog = fopen(mhFName.c_str(), "a"); + fprintf(fidlog,"\n"); + fprintf(fidlog,"%% Mh%dBlck%d ( %s %s )\n", (int) NewFile(b-1),b , __DATE__ , __TIME__ ); + fprintf(fidlog," \n"); + fprintf(fidlog," Number of simulations.: %d \n", currInitSizeArray);// (length(logpo2)) '); + fprintf(fidlog," Acceptation rate......: %f \n", jsux/currInitSizeArray); + fprintf(fidlog," Posterior mean........:\n"); + for (size_t i=0; i zeta_fwrd, zeta_back, zeta_mixed, zeta_static; + const mxArray *lli_mx = mxGetField(M_, 0, "lead_lag_incidence"); + MatrixConstView lli(mxGetPr(lli_mx), mxGetM(lli_mx), mxGetN(lli_mx), mxGetM(lli_mx)); + if (lli.getRows() != 3 || lli.getCols() != n_endo) + mexErrMsgTxt("Incorrect lead/lag incidence matrix"); + for (size_t i = 0; i < n_endo; i++) + { + if (lli(0, i) == 0 && lli(2, i) == 0) + zeta_static.push_back(i); + else if (lli(0, i) != 0 && lli(2, i) == 0) + zeta_back.push_back(i); + else if (lli(0, i) == 0 && lli(2, i) != 0) + zeta_fwrd.push_back(i); + else + zeta_mixed.push_back(i); + } + + double qz_criterium = *mxGetPr(mxGetField(options_, 0, "qz_criterium")); + double lyapunov_tol = *mxGetPr(mxGetField(options_, 0, "lyapunov_complex_threshold")); + double riccati_tol = *mxGetPr(mxGetField(options_, 0, "riccati_tol")); + size_t presample = (size_t) *mxGetPr(mxGetField(options_, 0, "presample")); + size_t console_mode = (size_t) *mxGetPr(mxGetField(options_, 0, "console_mode")); + + std::vector varobs; + const mxArray *varobs_mx = mxGetField(options_, 0, "varobs_id"); + if (mxGetM(varobs_mx) != 1) + mexErrMsgTxt("options_.varobs_id must be a row vector"); + size_t n_varobs = mxGetN(varobs_mx); + std::transform(mxGetPr(varobs_mx), mxGetPr(varobs_mx) + n_varobs, back_inserter(varobs), + std::bind2nd(std::minus(), 1)); + + if (data.getRows() != n_varobs) + mexErrMsgTxt("Data has not as many rows as there are observed variables"); + + std::vector estSubsamples; + estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1)); + + std::vector estParamsInfo; + fillEstParamsInfo(mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD, + estParamsInfo); + fillEstParamsInfo(mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD, + estParamsInfo); + fillEstParamsInfo(mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr, + estParamsInfo); + fillEstParamsInfo(mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr, + estParamsInfo); + fillEstParamsInfo(mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar, + estParamsInfo); + EstimatedParametersDescription epd(estSubsamples, estParamsInfo); + + // Allocate LogPosteriorDensity object + int info; + LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static, + qz_criterium, varobs, riccati_tol, lyapunov_tol, info); + + // Construct arguments of compute() method + Matrix steadyState(n_endo, 1); + mat::get_col(steadyState, 0) = VectorConstView(mxGetPr(mxGetField(oo_, 0, "steady_state")), n_endo, 1); + + Vector estParams2(n_estParams); + estParams2 = estParams; + Vector deepParams(n_param); + deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1); + Matrix Q(n_exo); + Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo); + + Matrix H(n_varobs); + const mxArray *H_mx = mxGetField(M_, 0, "H"); + if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0) + H.setAll(0.0); + else + H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs); + + // Construct MHMCMC Sampler + RandomWalkMetropolisHastings rwmh(estParams2.getSize()); + // Construct GaussianPrior drawDistribution m=0, sd=1 + GaussianPrior drawGaussDist01(0.0, 1.0, -INFINITY, INFINITY, 0.0, 1.0); + // get Jscale = diag(bayestopt_.jscale); + const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_"); + Matrix Jscale(n_estParams); + Matrix Dscale(n_estParams); + //Vector vJscale(n_estParams); + Jscale.setAll(0.0); + VectorConstView vJscale(mxGetPr(mxGetField(bayestopt_, 0, "jscale")), n_estParams, 1); + for (size_t i = 0; i < n_estParams; i++) + Jscale(i, i) = vJscale(i); + blas::gemm("N", "N", 1.0, D, Jscale, 0.0, Dscale); + +// Matrix mh_bounds(n_estParams,2); + + // Compute the MHMCMC loop draws + // and get get last line run in the last MH block sub-array + size_t lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams2, deepParams, data, Q, H, presample, info, + nMHruns, fblock, nBlocks, Dscale, drawGaussDist01, epd, resultsFileStem, console_mode); + + // Cleanups + for (std::vector::iterator it = estParamsInfo.begin(); + it != estParamsInfo.end(); it++) + delete it->prior; + + return lastMHblockArrayLine; + +} + +void +mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) +{ + if (nrhs != 7) + mexErrMsgTxt("logposterior: exactly seven arguments are required."); + if (nlhs != 1) + mexErrMsgTxt("logposterior: exactly one return argument is required."); + + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); +// Check and retrieve the arguments + + if (!mxIsDouble(prhs[0]) || mxGetN(prhs[0]) != 1) + mexErrMsgTxt("logposterior: First argument must be a column vector of double-precision numbers"); + + VectorConstView estParams(mxGetPr(prhs[0]), mxGetM(prhs[0]), 1); + + if (!mxIsDouble(prhs[1])) + mexErrMsgTxt("logposterior: Second argument must be a matrix of double-precision numbers"); + + MatrixConstView data(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]), mxGetM(prhs[1])); + + if (!mxIsChar(prhs[2])) + mexErrMsgTxt("logposterior: Third argument must be a character string"); + + char *mexext_mx = mxArrayToString(prhs[2]); + std::string mexext(mexext_mx); + mxFree(mexext_mx); + + size_t fblock = (size_t) mxGetScalar(prhs[3]); + size_t nBlocks = (size_t) mxGetScalar(prhs[4]); + VectorConstView nMHruns(mxGetPr(prhs[5]), mxGetM(prhs[5]), 1); + assert(nMHruns.getSize() == nBlocks); + + MatrixConstView D(mxGetPr(prhs[6]), mxGetM(prhs[6]), mxGetN(prhs[6]), mxGetM(prhs[6])); + + // Compute MCMC MH Draws and get last line run in the last MH block sub-array + size_t lastMHblockArrayLine = logMCMCposterior(estParams, data, mexext, fblock, nBlocks, nMHruns, D); + + *mxGetPr(plhs[0]) = (double) lastMHblockArrayLine; +}