markovchainmontecarlosimulationusingthemetropolisalgorithm-source.nb

simulation de montre carlo

(* Content-type: application/mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 6.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       145,          7]
NotebookDataLength[     24001,        631]
NotebookOptionsPosition[     22886,        590]
NotebookOutlinePosition[     23546,        615]
CellTagsIndexPosition[     23503,        612]
WindowFrame->Normal
ContainsDynamic->True *)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
Cell["\<\
Markov Chain Monte Carlo Simulation Using the Metropolis Algorithm\
\>", "Section",
 CellFrameColor->RGBColor[
  0.6449835965514611, 0.758632791638056, 0.2516823071641108],
 FontColor->RGBColor[
  0.6449835965514611, 0.758632791638056, 0.2516823071641108]],

Cell[BoxData[
 RowBox[{"Manipulate", "[", "\[IndentingNewLine]", 
  RowBox[{
   RowBox[{
    RowBox[{"SeedRandom", "[", "sr", "]"}], ";", "\[IndentingNewLine]", 
    RowBox[{"n", "=", "1000"}], ";", "\[IndentingNewLine]", 
    RowBox[{"x1", "=", 
     RowBox[{"-", "4"}]}], ";", "\[IndentingNewLine]", 
    RowBox[{"y1", "=", "9"}], ";", "\[IndentingNewLine]", 
    RowBox[{"xy", "=", 
     RowBox[{"Table", "[", 
      RowBox[{
       RowBox[{"{", 
        RowBox[{"0", ",", "0"}], "}"}], ",", 
       RowBox[{"{", 
        RowBox[{"i", ",", "1", ",", "n"}], "}"}]}], "]"}]}], ";", 
    "\[IndentingNewLine]", 
    RowBox[{
     RowBox[{"xy", "[", 
      RowBox[{"[", 
       RowBox[{"1", ",", "1"}], "]"}], "]"}], "=", "x1"}], ";", 
    "\[IndentingNewLine]", 
    RowBox[{
     RowBox[{"xy", "[", 
      RowBox[{"[", 
       RowBox[{"1", ",", "2"}], "]"}], "]"}], "=", "y1"}], ";", 
    "\[IndentingNewLine]", 
    RowBox[{"na", "=", "0"}], ";", "\[IndentingNewLine]", 
    RowBox[{"xypts", "=", 
     RowBox[{"Table", "[", 
      RowBox[{
       RowBox[{
        RowBox[{"xm", "=", 
         RowBox[{
          RowBox[{"xy", "[", 
           RowBox[{"[", 
            RowBox[{"All", ",", "1"}], "]"}], "]"}], "[", 
          RowBox[{"[", 
           RowBox[{"i", "-", "1"}], "]"}], "]"}]}], ";", 
        RowBox[{"ym", "=", 
         RowBox[{
          RowBox[{"xy", "[", 
           RowBox[{"[", 
            RowBox[{"All", ",", "2"}], "]"}], "]"}], "[", 
          RowBox[{"[", 
           RowBox[{"i", "-", "1"}], "]"}], "]"}]}], ";", 
        RowBox[{"ndistx", "=", 
         RowBox[{"NormalDistribution", "[", 
          RowBox[{"xm", ",", "proposalSigma"}], "]"}]}], ";", 
        RowBox[{"ndisty", "=", 
         RowBox[{"NormalDistribution", "[", 
          RowBox[{"ym", ",", "proposalSigma"}], "]"}]}], ";", 
        "\[IndentingNewLine]", 
        RowBox[{"ux", "=", 
         RowBox[{"RandomReal", "[", "ndistx", "]"}]}], ";", 
        "\[IndentingNewLine]", 
        RowBox[{"uy", "=", 
         RowBox[{"RandomReal", "[", "ndisty", "]"}]}], ";", 
        "\[IndentingNewLine]", 
        RowBox[{"p2", "=", 
         RowBox[{"pdf", "/.", 
          RowBox[{"{", 
           RowBox[{
            RowBox[{"x", "\[Rule]", "ux"}], ",", 
            RowBox[{"y", "\[Rule]", "uy"}]}], "}"}]}]}], ";", 
        RowBox[{"p1", "=", 
         RowBox[{"pdf", "/.", 
          RowBox[{"{", 
           RowBox[{
            RowBox[{"x", "\[Rule]", "xm"}], ",", 
            RowBox[{"y", "\[Rule]", "ym"}]}], "}"}]}]}], ";", 
        "\[IndentingNewLine]", 
        RowBox[{"r", "=", 
         RowBox[{"p2", "/", "p1"}]}], ";", 
        RowBox[{"u2", "=", 
         RowBox[{"RandomReal", "[", "]"}]}], ";", "\[IndentingNewLine]", 
        RowBox[{"If", "[", 
         RowBox[{
          RowBox[{"u2", "\[LessEqual]", "r"}], ",", 
          RowBox[{
           RowBox[{
            RowBox[{"xy", "[", 
             RowBox[{"[", "i", "]"}], "]"}], "=", 
            RowBox[{"{", 
             RowBox[{"ux", ",", "uy"}], "}"}]}], ";", 
           RowBox[{"na", "=", 
            RowBox[{"na", "+", "1"}]}]}], ",", 
          RowBox[{
           RowBox[{"xy", "[", 
            RowBox[{"[", "i", "]"}], "]"}], "=", 
           RowBox[{"{", 
            RowBox[{"xm", ",", "ym"}], "}"}]}]}], "]"}], ";", 
        "\[IndentingNewLine]", 
        RowBox[{"xy", "[", 
         RowBox[{"[", "i", "]"}], "]"}]}], ",", "\[IndentingNewLine]", 
       RowBox[{"{", 
        RowBox[{"i", ",", "2", ",", "n"}], "}"}]}], "]"}]}], ";", 
    RowBox[{"plgraphics", "=", 
     RowBox[{"Graphics", "[", 
      RowBox[{"Text", "[", 
       RowBox[{
        RowBox[{"Style", "[", 
         RowBox[{
          RowBox[{"\"\<proposal \[Sigma] = \>\"", "<>", 
           RowBox[{"ToString", "[", "proposalSigma", "]"}], "<>", 
           "\"\<\\nacceptance rate = \>\"", "<>", 
           RowBox[{"ToString", "[", 
            RowBox[{"100.", "*", 
             RowBox[{"na", "/", "n"}]}], "]"}], "<>", "\"\<%\>\""}], ",", 
          "16", ",", "Bold"}], "]"}], ",", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"-", "3"}], ",", "9"}], "}"}], ",", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"-", "1"}], ",", "0"}], "}"}], ",", 
        RowBox[{"{", 
         RowBox[{"1", ",", "0"}], "}"}]}], "]"}], "]"}]}], ";", 
    "\[IndentingNewLine]", 
    RowBox[{"points", "=", 
     RowBox[{"ListPlot", "[", 
      RowBox[{"xypts", ",", 
       RowBox[{"PlotRange", "\[Rule]", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{
            RowBox[{"-", "5"}], ",", "10"}], "}"}], ",", 
          RowBox[{"{", 
           RowBox[{
            RowBox[{"-", "5"}], ",", "10"}], "}"}]}], "}"}]}], ",", 
       RowBox[{"Frame", "\[Rule]", "True"}], ",", 
       RowBox[{"FrameLabel", "\[Rule]", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"Style", "[", 
           RowBox[{"\"\<x\>\"", ",", "Italic"}], "]"}], ",", 
          RowBox[{"Style", "[", 
           RowBox[{"\"\<y\>\"", ",", "Italic"}], "]"}]}], "}"}]}]}], "]"}]}], 
    ";", 
    RowBox[{"xpoints", "=", 
     RowBox[{"ListPlot", "[", 
      RowBox[{
       RowBox[{"xy", "[", 
        RowBox[{"[", 
         RowBox[{"All", ",", "1"}], "]"}], "]"}], ",", 
       RowBox[{"PlotRange", "\[Rule]", "All"}], ",", 
       RowBox[{"Frame", "\[Rule]", "True"}], ",", 
       RowBox[{"PlotLabel", "\[Rule]", 
        RowBox[{"Row", "[", 
         RowBox[{"{", 
          RowBox[{
           RowBox[{"Style", "[", 
            RowBox[{"\"\<x\>\"", ",", "Italic"}], "]"}], ",", 
           "\"\< trace\>\""}], "}"}], "]"}]}], ",", 
       RowBox[{"FrameLabel", "\[Rule]", 
        RowBox[{"{", 
         RowBox[{"\"\<iteration\>\"", ",", 
          RowBox[{"Style", "[", 
           RowBox[{"\"\<x\>\"", ",", "Italic"}], "]"}]}], "}"}]}], ",", 
       RowBox[{"ImagePadding", "\[Rule]", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{"35", ",", "10"}], "}"}], ",", 
          RowBox[{"{", 
           RowBox[{"35", ",", "25"}], "}"}]}], "}"}]}]}], "]"}]}], ";", 
    RowBox[{"ypoints", "=", 
     RowBox[{"ListPlot", "[", 
      RowBox[{
       RowBox[{"xy", "[", 
        RowBox[{"[", 
         RowBox[{"All", ",", "2"}], "]"}], "]"}], ",", 
       RowBox[{"PlotRange", "\[Rule]", "All"}], ",", 
       RowBox[{"Frame", "\[Rule]", "True"}], ",", 
       RowBox[{"PlotLabel", "\[Rule]", 
        RowBox[{"Row", "[", 
         RowBox[{"{", 
          RowBox[{
           RowBox[{"Style", "[", 
            RowBox[{"\"\<y\>\"", ",", "Italic"}], "]"}], ",", 
           "\"\< trace\>\""}], "}"}], "]"}]}], ",", 
       RowBox[{"FrameLabel", "\[Rule]", 
        RowBox[{"{", 
         RowBox[{"\"\<iteration\>\"", ",", 
          RowBox[{"Style", "[", 
           RowBox[{"\"\<y\>\"", ",", "Italic"}], "]"}]}], "}"}]}], ",", 
       RowBox[{"ImagePadding", "\[Rule]", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{"35", ",", "10"}], "}"}], ",", 
          RowBox[{"{", 
           RowBox[{"35", ",", "25"}], "}"}]}], "}"}]}]}], "]"}]}], ";", 
    "\[IndentingNewLine]", 
    RowBox[{"GraphicsGrid", "[", 
     RowBox[{"{", 
      RowBox[{
       RowBox[{"{", 
        RowBox[{"plgraphics", ",", 
         RowBox[{"Show", "[", 
          RowBox[{"plcont", ",", "points"}], "]"}]}], "}"}], ",", 
       RowBox[{"{", 
        RowBox[{"xpoints", ",", "ypoints"}], "}"}]}], "}"}], "]"}]}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"{", 
      RowBox[{"proposalSigma", ",", ".2", ",", "\"\<proposal \[Sigma]\>\""}], 
      "}"}], ",", ".2", ",", "10.2", ",", "1.0"}], "}"}], ",", 
   "\[IndentingNewLine]", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"{", 
      RowBox[{"sr", ",", "123", ",", "\"\<random seed\>\""}], "}"}], ",", "1",
      ",", "123456", ",", "1"}], "}"}], ",", 
   RowBox[{"TrackedSymbols", "\[RuleDelayed]", "Manipulate"}], ",", 
   "\[IndentingNewLine]", 
   RowBox[{"Initialization", "\[RuleDelayed]", "\[IndentingNewLine]", 
    RowBox[{"{", "\[IndentingNewLine]", 
     RowBox[{
      RowBox[{"pdf", "=", 
       RowBox[{
        RowBox[{"0.087", " ", 
         SuperscriptBox["\[ExponentialE]", 
          RowBox[{
           FractionBox["1", "2"], " ", 
           RowBox[{"(", 
            RowBox[{
             RowBox[{
              RowBox[{"-", 
               RowBox[{"(", 
                RowBox[{"x", "-", "4"}], ")"}]}], " ", 
              RowBox[{"(", 
               RowBox[{
                RowBox[{"0.595", 
                 RowBox[{"(", 
                  RowBox[{"x", "-", "4"}], ")"}]}], "-", 
                RowBox[{"0.238", " ", 
                 RowBox[{"(", 
                  RowBox[{"y", "-", "3"}], ")"}]}]}], ")"}]}], "-", 
             RowBox[{
              RowBox[{"(", 
               RowBox[{
                RowBox[{
                 RowBox[{"-", "0.238"}], 
                 RowBox[{"(", 
                  RowBox[{"x", "-", "4"}], ")"}]}], "+", 
                RowBox[{"0.595", " ", 
                 RowBox[{"(", 
                  RowBox[{"y", "-", "3"}], ")"}]}]}], ")"}], "  ", 
              RowBox[{"(", 
               RowBox[{"y", "-", "3"}], ")"}]}]}], ")"}]}]]}], "+", 
        FractionBox[
         SuperscriptBox["\[ExponentialE]", 
          RowBox[{
           FractionBox["1", "2"], " ", 
           RowBox[{"(", 
            RowBox[{
             RowBox[{"-", 
              SuperscriptBox["x", "2"]}], "-", 
             SuperscriptBox["y", "2"]}], ")"}]}]], 
         RowBox[{"2", " ", "\[Pi]"}]]}]}], ";", 
      RowBox[{"plcont", "=", 
       RowBox[{"ContourPlot", "[", 
        RowBox[{"pdf", ",", 
         RowBox[{"{", 
          RowBox[{"x", ",", 
           RowBox[{"-", "4"}], ",", "10"}], "}"}], ",", 
         RowBox[{"{", 
          RowBox[{"y", ",", 
           RowBox[{"-", "4"}], ",", "10"}], "}"}], ",", 
         RowBox[{"PlotRange", "\[Rule]", 
          RowBox[{"{", 
           RowBox[{
            RowBox[{"{", 
             RowBox[{
              RowBox[{"-", "5"}], ",", "10"}], "}"}], ",", 
            RowBox[{"{", 
             RowBox[{
              RowBox[{"-", "5"}], ",", "10"}], "}"}], ",", 
            RowBox[{"{", 
             RowBox[{"0", ",", "0.25"}], "}"}]}], "}"}]}], ",", 
         RowBox[{"Contours", "\[Rule]", "15"}], ",", 
         RowBox[{"ContourShading", "\[Rule]", "False"}], ",", 
         RowBox[{"FrameLabel", "\[Rule]", 
          RowBox[{"{", 
           RowBox[{
            RowBox[{"Style", "[", 
             RowBox[{"\"\<x\>\"", ",", "Italic"}], "]"}], ",", 
            RowBox[{"Style", "[",