// ComPhys File: Gaussian.java import java.applet.*; import java.awt.*; import java.awt.event.*; import comphys.*; public class Gaussian extends Applet implements ActionListener, ItemListener, Runnable { // define Gaussian probability function double A = 1; double sigma = 1; double xCenter = 0; boolean twoGaussians; double A2 = 0.7; double sigma2 = 0.5; double xCenter2 = 6; double p (double x) { double prob = A * Math.exp( - (x - xCenter) * (x - xCenter) / (2 * sigma * sigma) ); if (twoGaussians) prob += A2 * Math.exp( - (x - xCenter2) * (x - xCenter2) / (2 * sigma2 * sigma2) ); return prob; } // Metropolis algorithm double x0 = -5; double x = x0; double delta = 1; boolean MetropolisStep () { double xtrial = x + delta * (2 * Math.random() - 1); double w = p(xtrial) / p(x); if (Math.random() <= w) { x = xtrial; return true; } return false; } long naccept = 0; long ntotal = 0; int stepsToDiscard = 0; void Metropolis () { for (int step = 1; step < stepsToDiscard; step++) MetropolisStep(); if (MetropolisStep()) ++naccept; ++ntotal; } // histogram int N = 401; int[] histogram = new int[N]; int hMax; void initializeHistogram () { for (int n = 0; n < N; n++) histogram[n] = 0; hMax = 0; x = x0; } double xMin = -10; double xMax = 10; int pN = N * 5 / 6; int[] probability = new int[N]; void computeProbability () { double pMax = 0; for (int n = 0; n < N; n++) { double x = xMin + n * (xMax - xMin) / (double) N; double prob = p(x); if (p(x) > pMax) pMax = prob; } for (int n = 0; n < N; n++) { double x = xMin + n * (xMax - xMin) / (double) N; probability[n] = (int) (p(x) / pMax * pN); } } int xWalker; int deltaWalker; void locateWalker () { xWalker = (int) ( (x - xMin) / (xMax - xMin) * N ); deltaWalker = (int) (delta / (xMax - xMin) * N); } void recordStep () { locateWalker(); if (xWalker >= 0 && xWalker < N) if (++histogram[xWalker] > hMax) hMax = histogram[xWalker]; } class Picture extends Canvas { int pixels = N; Picture () { setSize(pixels, pixels); setBackground(new Color(255, 250, 240)); } public void paint (Graphics g) { update(g); } Image offScreen; Graphics osg; Image walker; public void update (Graphics g) { if (offScreen == null) { offScreen = createImage(pixels, pixels); osg = offScreen.getGraphics(); try { walker = getImage(getDocumentBase(), "walker.gif"); } catch (Exception e) { walker = Toolkit.getDefaultToolkit().getImage("walker.gif"); } } osg.clearRect(0, 0, pixels, pixels); double scale = Math.max(hMax, N) / (double) pN; int N0 = N - N / 12; osg.setColor(Color.red); if (xWalker >= 0 && xWalker < N) { int yWalker = N0 - probability[xWalker]; osg.fillOval(xWalker - 4, yWalker - 4, 8, 8); osg.drawImage(walker, xWalker - 15, yWalker - 20, this); } osg.setColor(Color.lightGray); for (int n = 0; n < N; n++) osg.drawLine(n, N0, n, N0 - (int) (histogram[n] / scale)); osg.setColor(Color.green); for (int n = 1; n < N; n++) osg.drawLine(n - 1, N0 - probability[n - 1], n, N0 - probability[n]); osg.setColor(Color.red); int dy = 16; osg.fillRect(xWalker - deltaWalker, N0 + dy, 2 * deltaWalker, 8); int dx = 16; dy = 12; int y = N0 + dy; int x = N / 2; osg.setColor(Color.blue); osg.drawLine(x, N0 - 2, x, N0 + 2); osg.drawString(" 0.0", x - dx, y); x = N / 4; osg.drawLine(x, N0 - 2, x, N0 + 2); osg.drawString("-5.0", x - dx, y); x = 3 * N / 4; osg.drawLine(x, N0 - 2, x, N0 + 2); osg.drawString("+5.0", x - dx, y); double accept = naccept; if (ntotal > 0) accept /= ntotal; osg.fillRect(10, 8, (int) (80 * accept), 8); osg.drawRect(10, 8, 80, 8); osg.drawString("" + (int) (accept * 100) + "% acceptance," + " Total steps = " + ntotal, 100, 16); g.drawImage(offScreen, 0, 0, null); } } Picture picture; TextField x0TextField; TextField deltaTextField; TextField discardTextField; TextField skipTextField; int skip; Checkbox twoGaussiansCheckbox; Button runButton; Button thermButton; Button resetButton; public void init () { initializeHistogram(); computeProbability(); locateWalker(); add(picture = new Picture()); Panel p = new Panel(); p.setLayout(new GridLayout(6, 2)); p.add(new Label("Start point x0 = ")); p.add(x0TextField = new TextField("" + x0)); x0TextField.addActionListener(this); p.add(new Label("Step size delta = ")); p.add(deltaTextField = new TextField("" + delta)); deltaTextField.addActionListener(this); p.add(new Label("Steps to discard = ")); p.add(discardTextField = new TextField("" + stepsToDiscard)); discardTextField.addActionListener(this); p.add(new Label("Skip plot steps = ")); p.add(skipTextField = new TextField("" + skip)); skipTextField.addActionListener(this); p.add(thermButton = new Button("Junk Therm Steps")); thermButton.addActionListener(this); p.add(twoGaussiansCheckbox = new Checkbox("Two peaks")); twoGaussiansCheckbox.addItemListener(this); p.add(runButton = new Button("Start")); runButton.addActionListener(this); p.add(resetButton = new Button("Reset")); resetButton.addActionListener(this); add(p); } Thread runThread; boolean running; public void actionPerformed (ActionEvent event) { if (event.getSource() == x0TextField) { try { x0 = new Double(event.getActionCommand()).doubleValue(); } catch (NumberFormatException nfe) { } x0TextField.setText("" + x0); x = x0; locateWalker(); } else if (event.getSource() == deltaTextField) { try { delta = new Double(event.getActionCommand()).doubleValue(); } catch (NumberFormatException nfe) { } deltaTextField.setText("" + delta); } else if (event.getSource() == discardTextField) { try { stepsToDiscard = Integer.parseInt(event.getActionCommand()); } catch (NumberFormatException nfe) { } discardTextField.setText("" + stepsToDiscard); } else if (event.getSource() == skipTextField) { try { skip = Integer.parseInt(event.getActionCommand()); } catch (NumberFormatException nfe) { } skipTextField.setText("" + skip); } else if (event.getSource() == thermButton) { naccept = ntotal = 0; initializeHistogram(); } else if (event.getSource() == runButton) { if (running) { running = false; runButton.setLabel("Start"); } else { running = true; runThread = new Thread(this); runThread.start(); runButton.setLabel("Stop"); } } else if (event.getSource() == resetButton) { if (running) { running = false; } naccept = ntotal = 0; initializeHistogram(); x = x0; locateWalker(); picture.repaint(); runButton.setLabel("Start"); } } public void itemStateChanged (ItemEvent itemEvent) { if (itemEvent.getSource() == twoGaussiansCheckbox) { twoGaussians = twoGaussiansCheckbox.getState(); computeProbability(); } } public void run () { runThread.setPriority(Thread.MIN_PRIORITY); while (running) { Metropolis(); recordStep(); picture.repaint(); for (int s = 0; s < skip; s++) { Metropolis(); recordStep(); } try { runThread.sleep(10); } catch (InterruptedException ie) { } } runThread = null; } public static void main (String[] args) { Gaussian gaussian = new Gaussian(); CPFrame aFrame = new CPFrame("The Gaussian distribution"); aFrame.add(gaussian); gaussian.init(); aFrame.setSize(700, 450); aFrame.setLocation(50, 50); aFrame.setVisible(true); } }