use of cern.jet.random.Gamma in project Gemma by PavlidisLab.
the class ComBat method plot.
/**
* Make diagnostic plots.
* FIXME: As in the original ComBat, this only graphs the first batch's statistics. In principle we can (and perhaps
* should) examine these plots for all the batches.
*
* @param filePrefix file prefix
*/
public void plot(String filePrefix) {
if (this.gammaHat == null)
throw new IllegalArgumentException("You must call 'run' first");
/*
* View the distribution of gammaHat, which we assume will have a normal distribution
*/
DoubleMatrix1D ghr = gammaHat.viewRow(0);
int NUM_HIST_BINS = 100;
Histogram gammaHatHist = new Histogram("GammaHat", NUM_HIST_BINS, ghr);
XYSeries ghplot = gammaHatHist.plot();
Normal rn = new Normal(this.gammaBar.get(0), Math.sqrt(this.t2.get(0)), new MersenneTwister());
Histogram ghtheoryT = new Histogram("Gamma", NUM_HIST_BINS, gammaHatHist.min(), gammaHatHist.max());
for (int i = 0; i < 10000; i++) {
double n = rn.nextDouble();
ghtheoryT.fill(n);
}
XYSeries ghtheory = ghtheoryT.plot();
File tmpfile;
try {
tmpfile = File.createTempFile(filePrefix + ".gammahat.histogram.", ".png");
ComBat.log.info(tmpfile);
} catch (IOException e) {
throw new RuntimeException(e);
}
try (OutputStream os = new FileOutputStream(tmpfile)) {
this.writePlot(os, ghplot, ghtheory);
/*
* View the distribution of deltaHat, which we assume has an inverse gamma distribution
*/
DoubleMatrix1D dhr = deltaHat.viewRow(0);
Histogram deltaHatHist = new Histogram("DeltaHat", NUM_HIST_BINS, dhr);
XYSeries dhplot = deltaHatHist.plot();
Gamma g = new Gamma(aPrior.get(0), bPrior.get(0), new MersenneTwister());
Histogram deltaHatT = new Histogram("Delta", NUM_HIST_BINS, deltaHatHist.min(), deltaHatHist.max());
for (int i = 0; i < 10000; i++) {
double invg = 1.0 / g.nextDouble();
deltaHatT.fill(invg);
}
XYSeries dhtheory = deltaHatT.plot();
tmpfile = File.createTempFile(filePrefix + ".deltahat.histogram.", ".png");
ComBat.log.info(tmpfile);
try (OutputStream os2 = new FileOutputStream(tmpfile)) {
this.writePlot(os2, dhplot, dhtheory);
}
} catch (IOException e) {
throw new RuntimeException(e);
}
}
Aggregations