Class for computing the Hessian of the negative free energy. More...
Public Member Functions | |
| FreeEnergyHessian (const DataForMCMC &d) | |
| Matrix< 4 > | hessian (const vector< Vector< 4 > > &spots, int spot) const |
| MT19937 & | get_rng () const |
Protected Attributes | |
| const vector< ImageRef > & | pixels |
| const vector< vector< double > > & | pixel_intensities |
| const double | mu_brightness |
| const double | sigma_brightness |
| const double | mu_blur |
| const double | sigma_blur |
| const double | variance |
| const int | samples |
| const int | sample_iterations |
| const Matrix< 3 > | A |
| const Vector< 3 > | pi |
| MT19937 & | rng |
Detailed Description
Class for computing the Hessian of the negative free energy.
Definition at line 608 of file multispot5.cc.
Constructor & Destructor Documentation
| FreeEnergyHessian::FreeEnergyHessian | ( | const DataForMCMC & | d | ) | [inline] |
Constructor.
- Parameters:
-
d All data required
Definition at line 614 of file multispot5.cc.
:DataForMCMC(d) { }
Member Function Documentation
| Matrix<4> FreeEnergyHessian::hessian | ( | const vector< Vector< 4 > > & | spots, |
| int | spot | ||
| ) | const [inline] |
Compute the Hessian.
- Parameters:
-
spots All spot positions spot spot to compute hessian for
Definition at line 622 of file multispot5.cc.
References SampledMultispot::compute_spot_intensity(), SampledMultispot::compute_spot_intensity_hessian(), diff_log_log_normal(), hess_log_log_normal(), SampledMultispot::GibbsSampler::next(), SampledMultispot::GibbsSampler::sample(), and SampledMultispot::GibbsSampler::sample_intensities().
Referenced by FitSpots::optimize_each_spot_in_turn_for_several_passes().
{
cout << "----Computing pure MCMC hessian\n";
const unsigned int nspots = spots.size();
const unsigned int nframes = pixel_intensities.size();
const unsigned int npixels = pixels.size();
cout << spot << " " << nspots << " " << nframes << " " << npixels << endl;
vector<vector<double> > spot_intensity; //[spot][pixel]
for(unsigned int i=0; i < nspots; i++)
spot_intensity.push_back(compute_spot_intensity(pixels, spots[i]));
vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(pixels, spots[spot]);
GibbsSampler sampler(pixel_intensities, spot_intensity, spots, A, pi, variance, sample_iterations);
//Compute derivative of log probability, summed (ie integrated)
Matrix<4> sum_hess1 = Zeros(spots.size());
Matrix<4> sum_hess2 = Zeros(spots.size());
Vector<4> sum_deriv = Zeros(spots.size());
for(int sample=0; sample < samples; sample++)
{
sampler.next(rng);
//Compute d log P(data | x, model) / d model, for a given sample
//And the hessian
Matrix<4> hess = Zeros(spots.size());
Vector<4> deriv = Zeros(spots.size());
for(unsigned int frame=0; frame < nframes; frame++)
{
for(unsigned int pixel=0; pixel < npixels; pixel++)
{
double e = pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel];
//Build up the derivative
if(sampler.sample()[spot][frame] == 0)
{
hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row();
deriv += e * get<1>(spot_hess_etc[pixel]);
}
}
}
hess[0][0] += hess_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness);
hess[1][1] += hess_log_log_normal(spots[spot][1], mu_blur, sigma_blur);
sum_hess1 += hess;
deriv[0] += diff_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness);
deriv[1] += diff_log_log_normal(spots[spot][1], mu_blur, sigma_blur);
sum_hess2 += deriv.as_col() * deriv.as_row();
sum_deriv = sum_deriv + deriv;
}
sum_hess1 /= (samples * variance);
sum_hess2 /= (samples * variance);
sum_deriv /= (samples * variance);
cout << sum_hess1 << endl;
cout << sum_hess2 << endl;
cout << sum_deriv.as_col() * sum_deriv.as_row() << endl;
cout << "......." << sum_deriv << endl;
//Put in the prior
//The derivative prior parts cancel out.
//Rather sensibly this means that the second derivatives can be
//computed without reference to the prior, and then the prior can
//be added in later, i.e.: P(data, parameters) = P(data|parameter) P(parameters)
//The second derivatives have been constructed to be diagonal
DiagonalMatrix<4> hess_prior = Zeros(spots.size());
cout << "sum of parts = \n" << sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row() << endl;
//TooN cannot currently add DiagonalMatrix to Matrix!!
//sum_hess.diagonal_slice() += hess_prior.diagonal_slice();
cout << "++++Done Computing pure MCMC hessian\n";
return sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row();
}
| MT19937& DataForMCMC::get_rng | ( | ) | const [inline, inherited] |
Definition at line 189 of file multispot5.cc.
Referenced by NegativeFreeEnergy::compute_with_mask(), and NegativeFreeEnergy::operator()().
{
return rng;
}
Member Data Documentation
const vector<ImageRef>& DataForMCMC::pixels [protected, inherited] |
Definition at line 179 of file multispot5.cc.
const vector<vector<double> >& DataForMCMC::pixel_intensities [protected, inherited] |
Definition at line 180 of file multispot5.cc.
const double DataForMCMC::mu_brightness [protected, inherited] |
Definition at line 181 of file multispot5.cc.
const double DataForMCMC::sigma_brightness [protected, inherited] |
Definition at line 181 of file multispot5.cc.
const double DataForMCMC::mu_blur [protected, inherited] |
Definition at line 181 of file multispot5.cc.
const double DataForMCMC::sigma_blur [protected, inherited] |
Definition at line 181 of file multispot5.cc.
const double DataForMCMC::variance [protected, inherited] |
Definition at line 182 of file multispot5.cc.
const int DataForMCMC::samples [protected, inherited] |
Definition at line 183 of file multispot5.cc.
const int DataForMCMC::sample_iterations [protected, inherited] |
Definition at line 183 of file multispot5.cc.
const Matrix<3> DataForMCMC::A [protected, inherited] |
Definition at line 184 of file multispot5.cc.
const Vector<3> DataForMCMC::pi [protected, inherited] |
Definition at line 185 of file multispot5.cc.
MT19937& DataForMCMC::rng [protected, inherited] |
Definition at line 186 of file multispot5.cc.
The documentation for this class was generated from the following file:
