 |
MRS
1.0
A C++ Class Library for Statistical Set Processing
|
|
|
MRS
1.0
A C++ Class Library for Statistical Set Processing
|
|
The Moore rejection sampler class for trans-dimensional targets over labeled metric spaces. More...
Public Member Functions | |
| MRSampler (Fobj &f, int max_n_boxes, double Alb, unsigned seed=gsl_rng_default_seed, bool use_f_scale=true) | |
| Initialised constructor. | |
| ~MRSampler () | |
| Destructor. | |
| double | getPALB () |
| Return lower bound on the acceptance prob. | |
| void | Refine (double Alb) |
| Do further partitioning until acc. prob. lower bound > Alb, then setup pdf. | |
| void | RefineUntil (unsigned int Desired_N_boxes) |
| Refine partition by bisections until Desired_N_boxes many boxes is reached, then setup pdf. | |
| void | Refine (int Nbisect) |
| Refine partition by doing Nbisect many bisections. | |
| int | get_n_boxes () |
| Return the number of boxes in RangeDomainSet. | |
| void | Print_Domain_Partition (std::ostream &out) |
| Print labeled boxes in domain partition DomainParts [C-XSC output format]. | |
| void | Output_Domain_Partition (std::ostream &out) |
| Print labeled boxes in domain partition DomainParts [naive TAB-delimited numeric only format]. | |
| std::ostream & | MRSoutput (std::ostream &os, const double eps=0) const |
| Print the RangeDomainSet in tab-delimited numeric only format. | |
| LabPnt | RejectionSampleOnce (int &tries) |
| Return one sample of labeled point via rejection sampling, if possible. | |
| void | RejectionSampleMany (int nRS, RSSample &theSample) |
| Draw nRS many samples of labeled points via rejection sampling, if possible. | |
| void | ImportanceSampleManyQuasi (int NSamples, bool residual, ISSample &theSample) |
| Importance sampling with Quasi random numbers -- [Ignore: experimental]. | |
| void | ImportanceSampleMany (int NSamples, bool residual, bool pseudoRNG, ISSample &theSample) |
| Importance sampling with Pseudo/Quasi random numbers. | |
| void | ImportanceSampleManyPseudo (int NSamples, bool residual, ISSample &theSample) |
| Importance sampling with Pseudo random numbers. | |
| void | PrintBoxes (int Nprint) |
| double | getIU () |
| double | getIL () |
| double | get_unscaled_IU () |
| double | getIUminusL () |
| double | getUmax () |
| double | getPAest () |
| double | get_wsum () |
| int | get_nsum () |
| double | get_wmax () |
| double | get_wmin () |
| void | updateIntegral () |
| update Integral using present partition | |
| void | updateUmax () |
| Updating the upper and lower bounds for the global maximum of target F. | |
| int | get_n_topologies () |
| int | get_nonresidual_samples () |
The Moore rejection sampler class for trans-dimensional targets over labeled metric spaces.
| std::ostream & MRSampler::MRSoutput | ( | std::ostream & | os, |
| const double | eps = 0 |
||
| ) | const |
Print the RangeDomainSet in tab-delimited numeric only format.
{
// do nothing if there is nothing in the set
if(!RangeDomainSet.empty())
{
RangedLabBoxSet::const_iterator it;
it = RangeDomainSet.begin();
double vol = Volume((it->LBox).Box);
while (vol>eps && it != RangeDomainSet.end())
{ // pull em out from top of pq
RangedLabBox theBox = *it;
ivector x = it->LBox.Box;
// label
os << (it->LBox.L);
// range enclosure
os << "\t" << Inf(it->BoxRE) << "\t" << Sup(it->BoxRE);
//then the box
for (int i = Lb(x); i <= Ub(x) ; i++)
{
os << "\t" << Inf(x[i]) << "\t" << Sup(x[i]);
}
os<<endl;
vol = Volume((it->LBox).Box);
it++;
} // end iteration through the set
} // end if set not empty
return os;
}
| void MRSampler::PrintBoxes | ( | int | Nprint | ) |
Print boxes with MATLAB.
{
Nprint = 0; //output through cout
int NBoxesToPrint=1000;
RangedLabBox theBox;
RangedLabBoxSet::const_iterator it = RangeDomainSet.begin ();
for (int i = 0; it != RangeDomainSet.end () && i < NBoxesToPrint; it++, i++)
{ // pull em out from top of pq
theBox = *it;
ivector x = it->LBox.Box;
//cout << "i: " << i << " Box: " << it->LBox.Box << " RE: :"
// << it->BoxRE << " Vol: " << it->BoxVol;
//cout << " IntDiam: " << it->BoxIntegral (UsingLogDensity) << endl;
cout << "rectangle('Position',["
<< Inf(x[1]) << "," << Inf(it->BoxRE) << "," << Sup(x[1])-Inf(x[1])
<< "," << Sup(it->BoxRE)-Inf(it->BoxRE)
<< "], 'FaceColor','b')" << endl;
} cout << endl;
}
| void MRSampler::RejectionSampleMany | ( | int | nRS, |
| RSSample & | theSample | ||
| ) |
Draw nRS many samples of labeled points via rejection sampling, if possible.
Draw nRS many rejection samples, if possible,.
RejectionSampleMany stores the samples in the RSSample object theSample.
{
LabPnt proposed_LPnt;
int total_tries = 0;
//ofstream psampout("MRS_proposed.samples");
for (; theSample.Samples.size () < unsigned(nRS);)
{
int proposed_index = static_cast<int>(gsl_ran_discrete
(rgsl, gslpdfstruct));
LabBox LBox = DomainParts[proposed_index];
rvector proposed_point =
F.DrawFromBoxPrior(LBox, DrawUnifUnitBox(rgsl, VecLen(LBox.Box)));
proposed_LPnt.Pnt = proposed_point;
proposed_LPnt.L = DomainParts[proposed_index].L;
//proposed_LPnt.Print(psampout);
real rand = gsl_rng_uniform (rgsl);
real Fprop;
if (rand > 1.0)
{
printf
("#proposed_index, UBox[proposed_index], height: %i %g %g \n",
proposed_index, _double (UBox[proposed_index]), _double (rand));
getchar ();
}
if (UsingLogDensity)
{
// real rany = rand * UBox[proposed_index];
// cout << "RSmany. rany: " << rany << " Ubox: "
// << UBox[proposed_index] << endl;
if (((RS_SQUEEZE) &&
(ln(rand) <= LoBox[proposed_index] - UBox[proposed_index]))
|| ln (rand) <= F(proposed_LPnt) - UBox[proposed_index])
{
// proposed_LPnt.Print(cout);
theSample.Samples.push_back (proposed_LPnt);
}
}
else
{
real rany = rand * UBox[proposed_index];
// less than lower bound, don't need to eval. function
if (((RS_SQUEEZE) && (rany <= LoBox[proposed_index]))
|| (rany <= (Fprop = F (proposed_LPnt))))
{
// proposed_LPnt.Print(cout);
theSample.Samples.push_back (proposed_LPnt);
}
}
// } // loop over tries
// if(tries == MAXNTRIES){ cout << "No RS accepted in "
// << MAXNTRIES << " tries. Refine partition?" << endl; }
total_tries++; // += tries;
// cout << "RS accepted, tries: " << theSample.Samples.size()
// << " " << total_tries << endl;
} // loop over samples
theSample.Nprop = total_tries;
theSample.Ntopologies = DomainLabelSet.size ();
theSample.LabelSet = DomainLabelSet;
theSample.n_Dim_Max = n_dim_max;
theSample.EnvelopeIntegral = getIU ();
// cerr << "In MRSample.SampleOnce. After " << MAXNTRIES
// << " proposals, none accepted. Acceptance prob. very low. " << endl;
// return proposed_LPnt;
}
| void MRSampler::updateUmax | ( | ) |
Updating the upper and lower bounds for the global maximum of target F.
For each box, first get upper bound U = Sup(range enclosure), and then Umax is the max over boxes of this and already have found during bisection Lmax, i.e. max of over boxes of Inf(range enclosure) Lmax, Umax are rigorous lower and upper bounds for global maximum of f each time U is eval'd for a box and found to be > current Umax, eval f at midpoint of box. Then max of these is fmid_max
{
RangedLabBox theBox;
real fmid_max = BIGNEGATIVE;
real f_scale_local;
bool f_scaleDone_local = false;
RangedLabBoxSet::const_iterator it = RangeDomainSet.begin ();
bool first = true;
for (; it != RangeDomainSet.end (); it++)
{
theBox = *it;
real U = Sup (theBox.BoxRE);
if (U > Umax || first)
{
Umax = U;
cout << "Umax: " << Umax << endl;
LabPnt lab_midpnt;
lab_midpnt.Pnt = mid (theBox.LBox.Box);
lab_midpnt.L = theBox.LBox.L;
real fmid = F (lab_midpnt);
if (fmid > fmid_max)
// f eval'd in middle of box with greatest U.
fmid_max = fmid;
first = false;
}
}
// Umax rigorous upper bound on f
// fMaxLB is lower bound on maximum of f
fMaxLB = (Lmax > fmid_max) ? Lmax : fmid_max;
if (UsingLogDensity)
{
// could use F - scale
f_scale_local = Umax - ln (UmaxMAX);
cout << "UmaxMAX, Umax, f_scale_local: " << UmaxMAX << " "
<< Umax << " " << f_scale_local << endl;
if (Umax - fMaxLB < LOGDIAMFMAX)
{
f_scaleDone_local = true;
}
}
else
{
// could use F/f_scale
f_scale_local = Umax / UmaxMAX;
if (Umax < exp (LOGDIAMFMAX) * fMaxLB)
{
f_scaleDone_local = true;
}
}
f_scale = (UseFScale)? f_scale_local: (UsingLogDensity)? 0: 1.0;
cout << "f_scale: " << f_scale << " " << f_scale_local << endl;
f_scaleDone = f_scaleDone_local;
cout << "bottom of updateUmax \n";
}