A simple example of coding and simulating a random process with the MCMC algorithm in JAVA, Part 1: The Gibbs Sampling

Suzuki Takashi (September 19, 2010 8:58 PM) | コメント(0)
<Caution>

･ This article uses three heavy java applets. So the images or the applets on your browser may act up under certaincircumstance.

･ How to deal with the problem is in waiting 30-60 seconds and not touching the keyboard and the mouse to display the correct images or applets on the browser.

This article is theoretically indebted to the scrupulous works of Yoshiyuki Wakui, Hiroshi Watanabe, and Takeshi Amemiya, which have a lot of lucid explanations and I'd very much appreciate their useful descriptions. In addition, I'd like to express my gratitude for some
helpful websites, such as http://codezine.jp/article/detail/116 as for how to make three-dimensional graph in JAVA and so on.

First I'd tentatively like to outline the situation we were put in with your agreement. A busy person generated normal random numbers with a computer. Then his phone rang. After he had answered the phone, someone
called out to him. After he had reluctantly wrapped up something urgent, he realized that he clean forgot the variance of the normal random variables. He barely remembered the mean and the form of the distribution. Unfortunately, he needed to generate the additional
normal random numbers which followed the almost same normal distribution.

1. The Data Set
The data set he generated is given below.

    n     1  2  3  4  5 
    y     0.58  -2.23  0.86  0.16  -0.20 
    n     6  7  8  9  10 
    y     0.54  1.52  -0.69  0.66  -1.75 
    n     11  12  13  14  15 
    y     -1.24  0.32  -0.58  -0.83  -0.74 
    n     16  17  18  19  20 
    y     -2.01  -0.71  -0.05  0.47  2.30 

The sample mean is about -0.18 and the unbiased sample variance is about 1.33. These random variables y are actually generated from the normal distribution N(0,1) and he knew the population mean of y was 0 but didn't know the population variance of σ2 was 1. That is to say we pretend we don't know the population variance of σ2 is 1 and keep carrying out the below experiment on the above condition.

2. The Estimation of The Prior Distribution
I supposed that the mean μ of the random variables followed the normal distribution, that the population mean of μ is 0, that the population variance of μ is in proportion to the population variance of the random
variables, and that the variance σ2 of the random variables followed the inverse gamma distribution. The example of the prior distribution is given below.
$$p\left( \theta \right) = p\left( \mu \right) p \left( \sigma^2 \right)$$

We substitute both p(μ) and p(σ2) in the above equation to find the prior distribution.
$$p\left( \mu \right) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp{\left( -\frac{\mu^2}{\frac{2\sigma^2}{m_0}} \right)} $$
$$p\left( \sigma^2 \right) = \frac{\beta^\alpha}{\Gamma\left(\alpha\right)} \left( \frac{1}{\sigma^2} \right)^{\alpha+1} \exp{ \left( -\frac{\beta}{\sigma^2} \right)} $$

We also substitute the parameters for both E(σ2) and Var(σ2) in the below equations to find the value of both α and β.
$$E\left( \sigma^2 \right) = \frac{\beta}{\alpha-1} \hspace{10mm}\alpha > 1 $$
$$Var\left( \sigma^2 \right) = \frac{\beta^2}{\left( \alpha+1 \right)^2 \left( \alpha-2 \right)} \hspace{10mm}\alpha > 2 $$

      source code: PriorDistribution.java

import java.awt.*;
import java.awt.event.*;
import java.applet.Applet;
import javax.swing.JLabel;
import javax.swing.JTextField;

public class PriorDistribution extends Applet implements ActionListener{

	// variable set
	private static final long serialVersionUID = 4049271078937921933L;
	final static int StartingX=50,StartingY=527;
	final static String[] rulerX = {"",""," 0.2",""," 0.6",""," 1.0",""," 1.4","","  1.8",""};
	final static String[] rulerY = {"",""," -1.0",""," -0.5","","   0",""," 0.5","","  1.0",""};
	final static String[] rulerZ = {"","","","",""," 0","","","","","","","",""};
	final int EndingX=StartingX+transformation(120,120,120).x;
	final int EndingY=StartingY+transformation(120,-120,-120).y;
	public int[] SmallestInnerPart=new int[EndingX+1];
	public int[] LargestInnerPart=new int[EndingX+1];
	public static double exp = 1.0;
	public static double var = 1.0;
	public static double ratio = 3.0;
	JTextField yx = new JTextField("1.0");
	JTextField yn = new JTextField("1.0");
	JTextField yr = new JTextField("3.0");
	JLabel label1 = new JLabel("The Prior Distribution", JLabel.CENTER);
	JLabel label2 = new JLabel("<html><body><font size=3>The Expectation of σ<sup>2</sup>(0<E(σ<sup>2</sup>)≤2.0)</font></body></html>", JLabel.CENTER);
	JLabel label3 = new JLabel("<html><body><font size=3>The Variance of σ<sup>2</sup> (0<Var(σ<sup>2</sup>)≤10<sup>8</sup>)</font></body></html>", JLabel.CENTER);
        JLabel label4 = new JLabel("<html><body><font size=3>    The Ratio of σ<sup>2</sup> to Var(μ) (0<r≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);

	public void init(){

		label1.setPreferredSize(new Dimension(400,24));
		label1.setFont(new Font("Serif",Font.BOLD,15));
		add(label1);
		label2.setPreferredSize(new Dimension(260,23));
		label2.setFont(new Font("Serif",Font.BOLD,11));
		add(label2);
		yx.setPreferredSize(new Dimension(120,23));
		add(yx);
		label3.setPreferredSize(new Dimension(260,23));
		label3.setFont(new Font("Serif",Font.BOLD,11));
		add(label3);
		yn.setPreferredSize(new Dimension(120,23));
		add(yn);
		label4.setPreferredSize(new Dimension(260,23));
		label4.setFont(new Font("Serif",Font.BOLD,11));
		add(label4);
		yr.setPreferredSize(new Dimension(120,23));
		add(yr);
		yx.addActionListener(this);
		yn.addActionListener(this);
		yr.addActionListener(this);
	}

	public void actionPerformed(ActionEvent e){
		if(e.getSource()==yx){
			exp = Double.valueOf(yx.getText()).doubleValue();
			if(exp <= 0 || exp > 2.0){exp=1.0;
			}
		}
		if(e.getSource()==yn){
			var = Double.valueOf(yn.getText()).doubleValue();
			if(var <= 0 || var > Math.pow(10.0,8)){var=1.0;
			}
		}
		if(e.getSource()==yr){
			ratio = Double.valueOf(yr.getText()).doubleValue();
			if(ratio <= 0 || ratio > Math.pow(10.0,8)){ratio=3.0;
			}
		}
		yx.setText(""+exp);
		yn.setText(""+var);
		yr.setText(""+ratio);
		repaint();
	}

	public void paint(Graphics g){

		Point dot,dot1,dot2,dot3,dot4,dot5,dot6;

		Graphics2D g2 = (Graphics2D)g;
		GradientPaint gp1 = new GradientPaint(0, 0, new Color(154,181,228), 0,470,new Color(225,232,245), true);
		g2.setPaint(gp1);
		g2.fillRect(0,0,460,447);
		super.paint(g);

		// x axis
		g.setColor(new Color(128,128,128));
		dot1=transformation(-120,-120,0);
                dot2=transformation(120,-120,0);
		g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
		for(int i=0;i<21;i++){
			dot1=transformation((120/10*i-120),-120,0);
			if(i % 4==0){
                     dot2=transformation((120/10*i-120),(120-10),0);
				g.setColor(new Color(128,128,128));
				g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
				g2.setFont(new Font ("SanSerif",Font.PLAIN,12));
				g.drawString(rulerX[i/2],dot2.x-48,dot2.y+15);
				if(i==8)
 					g.drawString("σ",dot2.x+28,dot2.y+30);
					g2.setFont(new Font ("SanSerif",Font.PLAIN,8));
				if(i==8) g.drawString("2",dot2.x+37,dot2.y+23);
			}
		}
         
		// y axis
		g.setColor(new Color(128,128,128));
		dot1=transformation(-120,-120,0);
		dot2=transformation(-120,120,0);
		dot3=transformation(-120,120,0);
		dot4=transformation(-120,120,210);
		dot5=transformation(120,-120,0);
		dot6=transformation(120,120,0);
		g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);

                g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
		g.drawLine(dot5.x,dot5.y,dot6.x,dot6.y);
		for(int i=0;i<210;i=i+1){
			dot3=transformation(-119,120,i);
			dot4=transformation(120,120,i);
			g.setColor(new Color(154,181,228));
                g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
		}
		g.setColor(new Color(128,128,128));
		for(int i=3;i<10;i++){
			dot3=transformation(-120,120,(12/4*(i+3)-15)*10);
			dot4=transformation(120,120,(12/4*(i+3)-15)*10);
                        g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
		}
		for(int i=0;i<21;i++){
			dot1=transformation(120,(12*i-120),0);
			if(i % 4==0){
				dot2=transformation((12+1)*10,(12*i-120),0);
				g.setColor(new Color(128,128,128));
				g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
				g.setFont(new Font ("SanSerif",Font.PLAIN,12));
				g.drawString(rulerY[i/2],dot2.x-10,dot2.y+13);
				if(i==8) g.drawString("μ",dot2.x+38,dot2.y+16);
			}
		}

		// z axis
		g.setColor(new Color(128,128,128));
		dot1=transformation(-120,-120,0);
		dot2=transformation(-120,-120,210);
		dot3=transformation(-120,119,210);
		dot4=transformation(-120,119,0);
		g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
		g.drawLine(dot2.x,dot2.y,dot3.x,dot3.y);
		g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
		for(int i=0;i<240;i=i+1){
			g.setColor(new Color(154,181,228));
			dot5=transformation(-120,(-119+i),210);
			dot6=transformation(-120,(-119+i),0);
			g.drawLine(dot5.x,dot5.y,dot6.x,dot6.y);
		}
		g.setColor(new Color(128,128,128));
		for(int i=2;i<10;i++){
			dot1=transformation(-120,-120,(30*(i+3)-150));
			dot2=transformation(-130,-120,(30*(i+3)-150));
			dot3=transformation(-120,-119,(30*(i+3)-150));
			dot4=transformation(-120,120,(30*(i+3)-150));
			g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
			g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
			g2.setFont(new Font ("SanSerif",Font.PLAIN,12));
			g.drawString(rulerZ[(i+3)],dot2.x-38,dot2.y+5);
		if(i==5) g.drawString("L",dot2.x-25,dot2.y-110);
		}

		// Coloring
		for(int i=0;i<=EndingX;i++){
			SmallestInnerPart[i]=EndingY;
			LargestInnerPart[i]=0;
		}
		dot=new Point();
		for(double x=120;x>=-120;x-=0.2)
			for(double y=-120;y<=120;y+=0.2){
				double transX=x/120+1;
				double transY=y/120;
				double transZ=function(transX,transY);
				double z=transZ*3.8;
				dot=transformation(x,y,z);
				if(dot.y<SmallestInnerPart[dot.x]){
					SmallestInnerPart[dot.x]=dot.y;
					g.setColor(gradation(z));
					g.drawRect(dot.x,dot.y,1,1);
				}
				if(dot.y>LargestInnerPart[dot.x]){
					LargestInnerPart[dot.x]=dot.y;
					g.setColor(new Color(128,128,128));
					g.drawRect(dot.x,dot.y,1,1);
				}
			}
		}

	public static double function(double x, double y){
		double z;
		double k0;
		double k1;
		double u;
		u=(Math.pow((2*Math.PI)*(x/ratio),(-0.5)))*Math.exp(-Math.pow((y-0),2)/(2*x/ratio));
	        k0=Math.pow(exp,2)/var+2;
		k1=Math.pow(exp,3)/var+exp;
		z=Math.pow(x,(-k0-1))*Math.exp((-k1)/(x))*u;
		if(x>=0){
			if(z>=0){
				if(u>=0)
		     			return 150*z;
				else return 0;
				}
			else return 0;
			}
		else	return 0;
	}

	public Color gradation(double z){
		int d,r,g,b;
		z=z*5.2;
   		if(z>=0) d=(int)z % 256;
		else     d=255-(-(int)z % 256);
		int c=(int)(d/85.334);
		switch(c){
			case 0:	r=79+76*d/86;
        				g=129+58*d/86;
		       		        b=189-100*d/86;
			break;
			case 1:	r=155+37*(d-86)/85;
					g=187-99*(d-86)/85;
					b=89-12*(d-86)/85;
			break;
			case 2:	r=192-113*(d-171)/84;
					g=88+41*(d-171)/84;
					b=77+112*(d-171)/84;
			break;
			default:	r=0; g=0; b=0;
			break;
		}
 		Color color=new Color(r,g,b);
		return color;
	}

	public Point transformation(double x,double y,double z){
		Point dot=new Point();
		dot.x=StartingX+(int)((x+120)+Math.cos(Math.toRadians(30))*(y+120)/2);
		dot.y=StartingY-(int)(Math.sin(Math.toRadians(30))*(y+120)/2+(z+120));
		return dot;
	}
      }

3. The Estimation of The Posterior Distribution
The example of the conditional likelihood function is
given below.
$$L\left( y|\theta \right) = \prod_{i=1}^{20}{ \frac{1}{\sqrt{2\pi\sigma^2}} \exp{\left( -\frac{\left( y_i-\mu\right)^2}{2\sigma^2} \right)} }$$

The posterior distribution is in proportion to the
product of the conditional likelihood function
multiplied by the prior distribution and its
example is given below.
$$p\left( \theta|y \right) \varpropto L\left( y|\theta \right) p \left( \theta \right) $$

      source code: PosteriorDistribution.java

import java.awt.*;
import java.awt.event.*;
import java.applet.Applet;
import javax.swing.JLabel;
import javax.swing.JTextField;

public class PosteriorDistribution extends Applet implements 
ActionListener{

    // variable set
    private static final long serialVersionUID = -6251476488227240006L;
    final static int StartingX=50,StartingY=527;
    final static String[] rulerX = {"",""," 0.2",""," 0.6",""," 1.0",""," 1.4","","  1.8",""};
    final static String[] rulerY = {"",""," -1.0",""," -0.5",""," 0",""," 0.5","","  1.0",""};
    final static String[] rulerZ = {"","","","",""," 0","","","","","","","",""};
    final int EndingX=StartingX+transformation(120,120,120).x;
    final int EndingY=StartingY+transformation(120,-120,-120).y;
    public int[] SmallestInnerPart=new int[EndingX+1];
    public int[] LargestInnerPart=new int[EndingX+1];
    public static double exp = 1.0;
    public static double var = 1.0;
    public static double ratio = 3.0;
    JTextField yx = new JTextField("1.0");
    JTextField yn = new JTextField("1.0");
    JTextField yr = new JTextField("3.0");
    JLabel label1 = new JLabel("The Posterior Distribution", JLabel.CENTER);
    JLabel label2 = new JLabel("<html><body><font size=3>The Expectation of σ<sup>2</sup> (0<E(σ<sup>2</sup>)≤2.0)</font></body></html>", variJLabel.CENTER);
    JLabel label3 = new JLabel("<html><body><font size=3>The Variance of σ<sup>2</sup> (0<Var(σ<sup>2</sup>)≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);
    JLabel label4 = new JLabel("<html><body><font size=3>The Ratio of σ<sup>2</sup> to Var(μ) (0<r≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);

    public void init(){

        label1.setPreferredSize(new Dimension(400,24));
        label1.setFont(new Font("Serif",Font.BOLD,15));
        add(label1);
        label2.setPreferredSize(new Dimension(260,23));
        label2.setFont(new Font("Serif",Font.BOLD,11));
        add(label2);
        yx.setPreferredSize(new Dimension(120,23));
        add(yx);
        label3.setPreferredSize(new Dimension(260,23));
        label3.setFont(new Font("Serif",Font.BOLD,11));
        add(label3);
        yn.setPreferredSize(new Dimension(120,23));
        add(yn);
        label4.setPreferredSize(new Dimension(260,23));
        label4.setFont(new Font("Serif",Font.BOLD,11));
        add(label4);
        yr.setPreferredSize(new Dimension(120,23));
        add(yr);
        yx.addActionListener(this);
        yn.addActionListener(this);
        yr.addActionListener(this);
    }

    public void actionPerformed(ActionEvent e){
        if(e.getSource()==yx){
            exp = Double.valueOf(yx.getText()).doubleValue();
            if(exp <= 0 || exp > 2.0){exp=1.0;
            }
                }
        if(e.getSource()==yn){
            var = Double.valueOf(yn.getText()).doubleValue();
            if(var <= 0 || var > Math.pow(10.0,8)){
                var=1.0;
            }
        }
        if(e.getSource()==yr){
            ratio = Double.valueOf(yr.getText()).doubleValue();
            if(ratio <= 0 || ratio >Math.pow(10.0,8)){ratio=3.0;
            }
        }
        yx.setText(""+exp);
        yn.setText(""+var);
        yr.setText(""+ratio);
        repaint();
    }

    public void paint(Graphics g){

        Point dot,dot1,dot2,dot3,dot4,dot5,dot6;

        Graphics2D g2 = (Graphics2D)g;
        GradientPaint gp1 = new GradientPaint(0, 0, new Color(154,181,228), 0,470,new Color(225,232,245), true);
        g2.setPaint(gp1);
        g2.fillRect(0,0,460,447);
        super.paint(g);

        // x axis
        g.setColor(new Color(128,128,128));
        dot1=transformation(-120,-120,0);
        dot2=transformation(120,-120,0);
        g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
        for(int i=0;i<21;i++){
            dot1=transformation((120/10*i-120),-120,0);
            if(i % 4==0){
                dot2=transformation((120/10*i-120),(-120-10),0);
                g.setColor(new Color(128,128,128));
                g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
                g2.setFont(new Font ("SanSerif",Font.PLAIN,12));
                g.drawString(rulerX[i/2],dot2.x-48,dot2.y+15);
                if(i==8)
                    g.drawString("σ",dot2.x+28,dot2.y+30);
                    g2.setFont(new Font("SanSerif",Font.PLAIN,8));
                if(i==8) g.drawString("2",dot2.x+37,dot2.y+23);
            }
        } 
 
        // y axis
        g.setColor(new Color(128,128,128));
        dot1=transformation(-120,-120,0);
        dot2=transformation(-120,120,0);
        dot3=transformation(-120,120,0);
        dot4=transformation(-120,120,210);
        dot5=transformation(120,-120,0);
        dot6=transformation(120,120,0);
        g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
        g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
        g.drawLine(dot5.x,dot5.y,dot6.x,dot6.y);
        for(int i=0;i<210;i=i+1){
            dot3=transformation(-119,120,i);
            dot4=transformation(120,120,i);
            g.setColor(new Color(154,181,228));
            g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
        }
        g.setColor(new Color(128,128,128));
        for(int i=3;i<10;i++){
            dot3=transformation(-120,120,(12/4*(i+3)-15)*10);
            dot4=transformation(120,120,(12/4*(i+3)-15)*10);
            g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
        }
        for(int i=0;i<21;i++){
            dot1=transformation(120,(12*i-120),0);
            if(i % 4==0){
                dot2=transformation((12+1)*10,(12*i-120),0);
                g.setColor(new Color(128,128,128));
                g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
                g.setFont(new Font ("SanSerif",Font.PLAIN,12));
                g.drawString(rulerY[i/2],dot2.x-10,dot2.y+13);
                if(i==8) g.drawString("μ",dot2.x+38,dot2.y+16);
           }
        }

        // z axis
        g.setColor(new Color(128,128,128));
        dot1=transformation(-120,-120,0);
        dot2=transformation(-120,-120,210);
        dot3=transformation(-120,119,210);
        dot4=transformation(-120,119,0);
        g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
        g.drawLine(dot2.x,dot2.y,dot3.x,dot3.y);
        g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
        for(int i=0;i<240;i=i+1){
            g.setColor(new Color(154,181,228));
            dot5=transformation(-120,(-119+i),210);
            dot6=transformation(-120,(-119+i),0);
            g.drawLine(dot5.x,dot5.y,dot6.x,dot6.y);
        }
        g.setColor(new Color(128,128,128));
        for(int i=2;i<10;i++){
            dot1=transformation(-120,-120,(30*(i+3)-150));
            dot2=transformation(-130,-120,(30*(i+3)-150));
            dot3=transformation(-120,-119,(30*(i+3)-150));
            dot4=transformation(-120,120,(30*(i+3)-150));
            g.drawLine(dot1.x,dot1.y,dot2.x,dot2.y);
            g.drawLine(dot3.x,dot3.y,dot4.x,dot4.y);
            g2.setFont(new Font ("SanSerif",Font.PLAIN,12));
            g.drawString(rulerZ[(i+3)],dot2.x-38,dot2.y+5);
        if(i==5) g.drawString("L",dot2.x-25,dot2.y-110);
        }

        // Coloring
        for(int i=0;i<=EndingX;i++){
            SmallestInnerPart[i]=EndingY;
            LargestInnerPart[i]=0;
        }
        dot=new Point();
        for(double x=120;x>=-120;x-=0.2)
            for(double y=-120;y<=120;y+=0.2){
                double transX=x/120+1;
                double transY=y/120;
                double transZ=function(transX,transY);
                double z=transZ*3.8;
                dot=transformation(x,y,z);
                if(dot.y<SmallestInnerPart[dot.x]){
                    SmallestInnerPart[dot.x]=dot.y;
                    g.setColor(gradation(z));
                    g.drawRect(dot.x,dot.y,1,1);
                }
                if(dot.y>LargestInnerPart[dot.x]){
                    LargestInnerPart[dot.x]=dot.y;
                    g.setColor(new Color(128,128,128));
                    g.drawRect(dot.x,dot.y,1,1);
                }
            }
        }

    public static double function(double x, double y){
        double z;
        double k0;
        double k1;
        double u;
        double d0=1;
        double[] d1= new double[21];
        double[] d2={0.58,-2.23,0.86,0.16,-0.20,0.54,1.52,-0.69,0.66,-1.75,-1.24,0.32,-0.58,-0.83,-0.74,-2.01,-0.71,-0.05,0.47,2.30};
        for(int i=0;i<20;i++){
            d1[i]=Math.pow(2*Math.PI*x,(-0.5))*Math.exp(-Math.pow((d2[i]-y),2)/(2*x));
            d0=d0*d1[i];
        }
        u=(Math.pow((2*Math.PI)*(x/ratio),(-0.5)))*Math.exp(-Math.pow((y-0),2)/(2*x/ratio));
        k0=Math.pow(exp,2)/var+2;
        k1=Math.pow(exp,3)/var+exp;
        z=160000000*Math.pow(x,(-k0-1))*Math.exp((-k1)/(x))*u;
        if(x>=0){
            if(z>=0){
                if(u>=0)
                         return 100000000*d0*z;
                else return 0;
                }
            else return 0;
            }
        else    return 0;
    }

    public Color gradation(double z){
        int d,r,g,b;
        z=z*5.2;
        if(z>=0) d=(int)z % 256;
        else     d=255-(-(int)z % 256);
        int c=(int)(d/85.334);
        switch(c){
            case 0:    r=79+76*d/86;
                          g=129+58*d/86;
                          b=189-100*d/86;
            break;
            case 1:    r=155+37*(d-86)/85;
                          g=187-99*(d-86)/85;
                          b=89-12*(d-86)/85;
            break;
            case 2:    r=192-113*(d-171)/84;
                          g=88+41*(d-171)/84;
                          b=77+112*(d-171)/84;
            break;
            default:    r=0; g=0; b=0;
            break;
        }
        Color color=new Color(r,g,b);
        return color;
    }

    public Point transformation(double x,double y,double z){
        Point dot=new Point();
        dot.x=StartingX+(int)((x+120)+Math.cos(Math.toRadians(30))*(y+120)/2);
        dot.y=StartingY-(int)(Math.sin(Math.toRadians(30))*(y+120)/2+(z+120));
        return dot;
    } 
}

4. The Gibbs Sampler Algorithm
I am most grateful to the helpful website, MY STRUGGLE, and the URL is given below.
http://vyshemirsky.blogspot.com/2007/11/sample-from-gamma-distribution-in-java.html
The example of the Gibbs sampler is given below.
$$\mu^*_{i+1}\sim N\left(\frac{n\bar{y}}{r+n},\frac{\sigma^{*^2}_i}{r}\right)$$
$$\frac{1}{\sigma^{*^2}_{i+1}}\sim Gamma\left(k,\theta_{i+1}\right)$$

We substitute the shape parameter k and the scale parameter θi+1 in the probability density function of
the inverse gamma distribution and I provisionally set σ1*2 to 1.
$$k=\frac{\left\{E\left(\sigma^2\right)\right\}^2}{Var\left(\sigma^2\right)}+\frac{5}{2}$$
$$\theta_{i+1}=\frac{2}{{\displaystyle2\left\{\frac{E\left(\sigma^2\right)^3}{Var\left(\sigma^2\right)}+E\left(\sigma^2\right)\right\}+\sum^{20}_{j=1}\left(y_j-\bar{y}\right)^2+\frac{rn}{r+n}\bar{y}^2+\left(r+n\right)\left(u^*_{i+1}-\frac{r\bar{y}}{r+n} \right)^2}}$$

      source code: GibbsSampler.java

import java.awt.*;
import java.awt.event.*;
import java.text.DecimalFormat;
import java.applet.Applet;
import javax.swing.JLabel;
import javax.swing.JTextField;

public class GibbsSampler extends Applet implements ActionListener{

	// variable set
	private static final long serialVersionUID = -729195055739028464L;
	public static double exp = 1.0;
	public static double var = 1.0;
	public static double m0 = 3.0;
	public static double u, d3, d4, d7, d8, d9, d10, d11, d12;
	public static double[] s1 = new double[10000];
	public static double[] s2 = new double[10000];
	public static double[] s3 = new double[1000];
	public static double[] s4 = new double[1000];
	public static double[] d1 = new double[10000];
	public static double[] d2 = new double[10000];
	public static double[] d5 = new double[10000];
	public static double[] d6 = new double[10000];
	JTextField yx = new JTextField("1.0");
	JTextField yn = new JTextField("1.0");
	JTextField yr = new JTextField("3.0");
	JLabel label1 = new JLabel("Gibbs Sampler Frequency Distribution",JLabel.CENTER);
	JLabel label2 = new JLabel("<html><body><font size=3>The Expectation of σ<sup>2</sup> (0<E(σ<sup>2</sup>)≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);
	JLabel label3 = new JLabel("<html><body><font size=3>The Variance of σ<sup>2</sup> (0<Var(σ<sup>2</sup>)≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);
	JLabel label4 = new JLabel("<html><body><font size=3>The Ratio of σ<sup>2</sup> to Var(μ) (0<r≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);

	public void init(){
		label1.setPreferredSize(new Dimension(416,24));
		label1.setFont(new Font("Serif",Font.BOLD,14));
		add(label1);
		label2.setPreferredSize(new Dimension(260,23));
		label2.setFont(new Font("Serif",Font.BOLD,11));
		add(label2);
		yx.setPreferredSize(new Dimension(80,23));
		add(yx);
		label3.setPreferredSize(new Dimension(260,23));
		label3.setFont(new Font("Serif",Font.BOLD,11));
		add(label3);
		yn.setPreferredSize(new Dimension(80,23));
		add(yn);
		label4.setPreferredSize(new Dimension(260,23));
		label4.setFont(new Font("Serif",Font.BOLD,11));
		add(label4);
		yr.setPreferredSize(new Dimension(80,23));
		add(yr);
		yx.addActionListener(this);
		yn.addActionListener(this);
		yr.addActionListener(this);
	}

	public void actionPerformed(ActionEvent e){
		if(e.getSource()==yx){
			exp=Double.valueOf(yx.getText()).doubleValue();
			if(exp <= 0 || exp > Math.pow(10.0,8)){
				exp = 1.0;
			}
		}
		if(e.getSource()==yn){
			var = Double.valueOf(yn.getText()).doubleValue();
			if(var <= 0 || var > Math.pow(10.0,8)){
				var = 1.0;
			}
		}
		if(e.getSource()==yr){m0 = Double.valueOf(yr.getText()).doubleValue();
			if(m0 <= 0 || m0 > Math.pow(10.0,8))
				{m0 = 3.0;}
		}
		yx.setText(""+exp);
		yn.setText(""+var);
		yr.setText(""+m0);
		repaint();
	}

	public void paint(Graphics g){
        double n = 20;
        double n0 = 2*Math.pow(exp,2)/var+4;
        double s0 = 2*(Math.pow(exp,3)/var+exp)/n0;
        double u0 = 0;
        double q = 25.31;
        double ExpX = -0.18;
        double m1 = m0+n;
        double n1 = n0+n;
        double n1s1 = n0*s0+q+(m0*n/(m0+n))*(ExpX-u0)*(ExpX-u0);
        double u1 = (n*ExpX+m0*u0)/(m0+n);
		s1[0]=1;
		s2[0]=normal.sampleNormal(u1,(1/s1[0])/m1);
		for(int i=0;i<9999;i++){
			s1[i+1]=gamma.sampleGamma((n1+1)/2,1/(n1s1/2+m1*Math.pow(s2[i]-u1,2)/2));
			s2[i+1]=normal.sampleNormal(u1,(1/s1[i+1])/m1);
		}

		d1[5000]=1/s1[5000];
		d2[5000]=s2[5000];
		for(int i=5001;i<10000;i++){
			d1[i]=d1[i-1]+1/s1[i];
			d2[i]=d2[i-1]+s2[i];
		}
		d3=d1[9999]/5000;
		d4=d2[9999]/5000;

		d5[4999]=0;
		d6[4999]=0;

		for(int i=5000;i<10000;i++){
			d5[i]=Math.pow((1/s1[i]-d3),2)+d5[i-1];
			d6[i]=Math.pow((s2[i]-d4),2)+d6[i-1];
		}

		d7=d5[9999]/(5000-1);
		d8=d6[9999]/(5000-1);

		// variable set
		int[] xxxx = new int[1000];
		int[] yyyy = new int[1000];
		int[] xxxx1 = new int[1000];
		int[] yyyy1 = new int[1000];
		double	minss1 = 1/s1[5000];
		double	maxss1 = 1/s1[5000];
		double	minss2 = s2[5000];
		double	maxss2 = s2[5000];

		for(int i=5000;i<10000;i++){
			if(1/s1[i]>maxss1){
				maxss1 = 1/s1[i];
			}
		}

		for(int i=5000;i<10000;i++){
			if(1/s1[i]<minss1){
				minss1 = 1/s1[i];
			}
		}

		for(int i=5000;i<10000;i++){
			if(s2[i]>maxss2){
				maxss2 = s2[i];
			}
		}

		for(int i=5000;i<10000;i++){
			if(s2[i]<minss2){
				minss2 = s2[i];
			}
		}
		for(int i=0;i<1000;i++){
			s3[i]=0;
			s4[i]=0;
		}

		double aaa=(maxss1-minss1)/1000;
		double bbb=(maxss2-minss2)/1000;
		for(int i=5000;i<10000;i++){
			for(int j=0;j<1000;j++){
				if(1/s1[i] >= (minss1+aaa*j) && 1/s1[i] < minss1+aaa*(j+1))
					s3[j]=(s3[j]+1.0);
			}
		}

		for(int i=5000;i<10000;i++){
			for(int j=0;j<1000;j++){
				if(minss2+bbb*j<=s2[i] && s2[i]<minss2+bbb*(j+1))
					s4[j]=s4[j]+1;
			}
		}

		// variable set
		double minss3 = s3[0];
		double maxss3 = s3[0];
		double minss4 = s4[0];
		double maxss4 = s4[0];

		for(int i=0;i<1000;i++){
			if(maxss3<s3[i]){
				maxss3 = s3[i];
			}
		}

		for(int i=0;i<1000;i++){
			if(s3[i]<minss3){
				minss3 = s3[i];
			}
		}

		for(int i=0;i<1000;i++){
			if(s4[i]>maxss4){
				maxss4 = s4[i];
			}
		}

		for(int i=0;i<1000;i++){
			if(s4[i]<minss4){
				minss4 = s4[i];
			}
		}

		double yyy = maxss3-minss3;
		double yyy1 = maxss4-minss4;

		for(int i=0;i<1000;i++){
			xxxx[i] = (int)(i*351/1000)+33;
			yyyy[i] = 474-(int)((s3[i]-minss3)*351/yyy);
			xxxx1[i] = (int)(i*351/1000)+33;
			yyyy1[i] = 474-(int)((s4[i]-minss4)*351/yyy1);
		}

		Graphics2D g2 = (Graphics2D)g;
		GradientPaint gp1 = new GradientPaint(0, 0, new Color(154,181,228), 0,470,new Color(225,232,245), true);
		g2.setPaint(gp1);
		g2.fillRect(0,0,416,503);
		super.paint(g);
		GradientPaint gp2 = new GradientPaint(0, 33, new Color(225,232,245), 0,351,new Color(154,181,228), true);
		g2.setPaint(gp2);
		g2.fillRect(33,123,354,351);
		GradientPaint gp3 = new GradientPaint(33, 227, new Color(79,129,189), 33,418,new Color(65,105,225), true);
		g2.setPaint(gp3);
		for (int i=0;i<1000;i++){
  			g2.drawLine(xxxx1[i], 474, xxxx1[i], yyyy1[i]);
		}
		GradientPaint(33, 227, new Color(192,88,77), 33,418,new Color(160,82,45), true);
		g2.setPaint(gp4);
		for (int i=0;i<1000;i++){
  			g2.drawLine(xxxx[i], 474, xxxx[i], yyyy[i]);
		}
		DecimalFormat df = new DecimalFormat("0.00");
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.setColor (new Color(0,0,0));
		g2.drawString("n = 5000",280,157);
		g2.setColor (new Color(0,0,0));
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("E(μ",280,177);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("*",300,176);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("|y) =",305,177);
		g2.drawString(df.format(d4),335,177);
		g2.drawString("Var(μ",280,197);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("*",312,196);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("|y) =",317,197);
		g2.drawString(df.format(d8),347,197);
		g2.setColor (new Color(0,0,0));
		g2.drawString("E(σ",280,217);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("*",300,216);
		g2.setFont(new Font ("Serif",Font.BOLD,8));
		g2.drawString("2",305,212);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("|y) =",310,217);
		g2.drawString(df.format(d3),340,217);
		g2.drawString("Var(σ",280,237);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("*",312,236);
		g2.setFont(new Font ("Serif",Font.BOLD,8));
		g2.drawString("2",317,232);
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.drawString("|y) =",322,237);
		g2.drawString(df.format(d7),352,237);
		g2.drawString("E(y) = -0.18",280,257);
		g2.drawString("Var(y) = 1.33",280,277);
		g2.setColor (new Color(0,0,0));
		g2.drawString("F",15,143);
		g2.setColor (new Color(79,129,189));
		g2.drawString("μ",345,488);
		g2.setColor (new Color(0,0,0));
		g2.drawString(",",353,488);
		g2.setColor (new Color(192,88,77));
		g2.drawString("σ",360,488);
		g2.setFont(new Font ("Serif",Font.BOLD,8));
		g2.drawString("2",371,483);
	}
}

         source code: normal.java

import java.util.*;

public class normal {
    private static Random rnd = new Random(Calendar.getInstance().
    getTimeInMillis()+Thread.currentThread().getId());
    private static double u;
    public static double sampleNormal(double e0, double e1) {
        double e2;
        e2=Math.pow(e1,(0.5));
        u=e2*rnd.nextGaussian()+e0;
        return u;
    }
}

            gamma.java

import java.util.*;

// See the below website, copy and paste the class "Samplers".
// Then modify the code as necessary.
//
// http://vyshemirsky.blogspot.com/2007/11/
// sample-from-gamma-distribution-in-java.html
//
// I express my appreciation for MY STRUGGLE.

In conclusion, the posterior distribution is roughly in such a position that the prior distribution moves from (μ, σ2) to (μ, σ2+0.5) because the posterior distribution is in proportion to the product of the conditional likelihood function and the prior distribution and the posterior distribution are the conjugate distributions. I also presume the assumption that I set the both expectation of σ2 and the variance of σ2 to small numbers reduces the value of both Var(μ*|y) and Var(σ*2|y).

Finally I am happy to assist you in coding and simulating a random process with the MCMC algorithm in JAVA.

カテゴリ:

コメントする

My Photo
プロフィール！
2016･11･15 改訂
spacer01
rssspacer01foaf
spacer01
atom.xml
spacer01

この記事について

このページは、Suzuki TakashiがSeptember 19, 2010 8:58 PMに書いた記事です。

ひとつ前の記事は「調査は続行ながらも、ひとまず区切りか。」です。

次の記事は「修正に次ぐ修正。」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。

June 2024

Sun Mon Tue Wed Thu Fri Sat
1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30

カテゴリ

月別 アーカイブ

OpenID対応しています OpenIDについて