Skip to content

Instantly share code, notes, and snippets.

@JEM-Mosig

JEM-Mosig/blog-mackay_gaussian_2006.md

Last active April 2, 2019 10:26
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save JEM-Mosig/39ac1a9f88d9bd7bc18c168aeb30dca1 to your computer and use it in GitHub Desktop.
Save JEM-Mosig/39ac1a9f88d9bd7bc18c168aeb30dca1 to your computer and use it in GitHub Desktop.
Download ZIP
blog-mackay_gaussian_2006
Raw
blog-mackay_gaussian_2006.md

Wolframscripts that generate figures in my blog post

The scripts in this gist generate the animated figures in my blog post mackay_gaussian_2006. You need Wolfram Mathematica 11 or higher to run these scripts. For example,

wolframscript <name>.wls

will generate the figure <name>.png in the current directory, where <name> is the name of the script.

Raw
mackay_gaussian_2006-bi-posterior.wls
#!/usr/bin/env wolframscript
(* ::Package:: *)
(* ::Text:: *)
(*This script generates the second animated figure on my blog post "mackay_gaussian_2006".*)
Print["https://jem-mosig.com/2019/03/mackay_gaussian_2006/"]
(* ::Text:: *)
(*Generate all frames of the animation, one sample from the bi-variate normal distribution for each frame.*)
Print["Generating figure. This may take a minute..."]
frames = BlockRandom[
SeedRandom[2019];
With[{x = {2}, \[Sigma] = 1, l = 1, x0 = 1, y0 = -1},
Module[{\[Mu], K, Kinv, A1, A2, B, \[Mu]post, Kpost},
\[Mu] = ConstantArray[0, Length[x] + 1];
K = N[Outer[\[Sigma]^2 Exp[-(#1 - #2)^2/(2 l^2)] &,
Join[{x0}, x], Join[{x0}, x]], 40];
Kinv = Inverse[K];
A1 = Kinv[[1 ;; 1, 1 ;; 1]];
B = Kinv[[1 ;; 1, 2 ;;]];
A2 = Kinv[[2 ;;, 2 ;;]];
Assert[
AllTrue[
Flatten@Chop[Kinv - ArrayFlatten[{{A1, B}, {Transpose[B], A2}}]], # == 0 &
]
];
\[Mu]post = -Inverse[A2].Transpose[B].{{y0}} // Flatten;
Kpost = Inverse[A2];
Quiet[
With[{
dist = MultinormalDistribution[\[Mu]post, Kpost],
ah =
Graphics[{Line[{{{0, -1}, {0, 1}}, {{-1, 0}, {0, 0}}}]}],
contourPlot = ContourPlot[
PDF[MultinormalDistribution[\[Mu], K], {y1, y2}],
{y1, -4, 4}, {y2, -4, 4},
Contours -> 10,
ScalingFunctions -> "Log",
ColorFunction -> GrayLevel,
ImageSize -> Small,
FrameLabel -> {"\!\(\*SubscriptBox[\(y\), \(1\)]\)",
"\!\(\*SubscriptBox[\(y\), \(2\)]\)"},
Epilog -> {Dashed, ColorData[112, 1],
InfiniteLine[{y0, 0}, {0, 1}]}
]
},
Table[
With[{y = RandomVariate[dist]},
Row[{
Show[{contourPlot,
Graphics[{ColorData[112, 2], PointSize[Large],
Point[{y0, y[[1]]}]}]}
],
ListPlot[
Prepend[Transpose[{x, y}], {x0, y0}],
Joined -> True, InterpolationOrder -> 1,
PlotStyle -> Directive[Thick, ColorData[112, 2]],
PlotRange -> {{0.8, 2.2}, {-3, 3}},
FrameTicks -> {{Automatic,
Automatic}, {{{1,
"\!\(\*SubscriptBox[\(y\), \(1\)]\)"}, {2,
"\!\(\*SubscriptBox[\(y\), \(2\)]\)"}}, None}},
Axes -> False, Frame -> True,
Epilog -> {
Black,
PointSize[Large],
Point[{x0, y0}]
},
Prolog -> {
Black, AbsoluteThickness[2],
PointSize[Large],
ColorData[112, 2],
Point[Transpose[{x, y}]],
Black,
Arrowheads[{{.025, 1., ah}, {-.025, 0., ah}}],
With[{i = 1},
Arrow[{{x[[i]], \[Mu]post[[i]] -
Sqrt@Kpost[[i, i]]}, {x[[i]], \[Mu]post[[i]] +
Sqrt@Kpost[[i, i]]}}]
],
Gray, Arrowheads[{{.025, 1., ah}, {-.025, 0., ah}}],
With[{i = 2},
Arrow[{{2 + 0.05, 0 - Sqrt@K[[i, i]]}, {2 + 0.05,
0 + Sqrt@K[[i, i]]}}]
]
},
ImageSize -> Medium
]
}]
],
100
]
],
{MultinormalDistribution::posdefprm, Transpose::nmtx}
]
]
]
];
(* ::Text:: *)
(*Sort the frames to appear visually more pleasing.*)
ordering = Quiet[
Last@FindShortestTour@Table[
FirstCase[frame[[1,1]],Point[pts_]:>pts,None,Infinity],
{frame,frames}
],
{InterpolatingFunction::dmval}
];
frames = frames[[ordering]];
Print@Export[
FileNameJoin[{"mackay_gaussian_2006-bi-posterior.gif"}],
frames,
AnimationRepetitions->Infinity,
"DisplayDurations"->.1
]
Raw
mackay_gaussian_2006-bi-prior.wls
#!/usr/bin/env wolframscript
(* ::Package:: *)
(* ::Text:: *)
(*This script generates the first animated figure on my blog post "mackay_gaussian_2006".*)
Print["https://jem-mosig.com/2019/03/mackay_gaussian_2006/"]
(* ::Text:: *)
(*Generate all frames of the animation, one sample from the bi-variate normal distribution for each frame.*)
Print["Generating figure. This may take a minute..."]
frames = BlockRandom[
SeedRandom[2019];
With[{x={1,2}, \[Sigma]=1, l=1},
Module[{\[Mu],K,Kinv,A1,A2,B,\[Mu]post,Kpost},
\[Mu] = ConstantArray[0,Length[x]];
K = N[Outer[\[Sigma]^2 Exp[-(#1-#2)^2/(2l^2)]&,x,x],40];
Kinv = Inverse[K];
Quiet[
With[{
dist = MultinormalDistribution[\[Mu],K],
ah = Graphics[{Line[{{{0,-1},{0,1}}}]}],
contourPlot = ContourPlot[
PDF[MultinormalDistribution[\[Mu],K],{y1,y2}],
{y1,-4,4},{y2,-4,4},
Contours->10,
ScalingFunctions->"Log",
ColorFunction->GrayLevel,
ImageSize->Small,
FrameLabel->{"\!\(\*SubscriptBox[\(y\), \(1\)]\)","\!\(\*SubscriptBox[\(y\), \(2\)]\)"}
]
},
Table[
With[{y=RandomVariate[dist]},
Row[{
Show[{contourPlot,Graphics[{ColorData[112,2],PointSize[Large],Point[y]}]}],
Spacer[10],
ListPlot[
Transpose[{x,y}],
Joined->True,
InterpolationOrder->1,
PlotStyle->Directive[Thick,ColorData[112,2]],
PlotRange->{{0.8,2.2},{-3,3}},
FrameTicks->{{Automatic,Automatic},{{{1,"\!\(\*SubscriptBox[\(y\), \(1\)]\)"},{2,"\!\(\*SubscriptBox[\(y\), \(2\)]\)"}},None}},
Axes->False,Frame->True,
Prolog->{
ColorData[112,2],PointSize[Large],Point[Transpose[{x,y}]],
Black,AbsoluteThickness[2],
Arrowheads[{{.025,1.,ah},{-.025,0.,ah}}],
Table[Arrow[{{x[[i]],\[Mu][[i]]-Sqrt@K[[i,i]]},{x[[i]],\[Mu][[i]]+Sqrt@K[[i,i]]}}],{i,2}]
},
ImageSize->Medium
]
}]
],
100
]
],
{MultinormalDistribution::posdefprm,Transpose::nmtx}
]
]
]
];
(* ::Text:: *)
(*Sort the frames to appear visually more pleasing.*)
ordering = Quiet[
Last@FindShortestTour@Table[
FirstCase[frame[[1,1]],Point[pts_]:>pts,None,Infinity],
{frame,frames}
],
{InterpolatingFunction::dmval}
];
frames = frames[[ordering]];
(* ::Text:: *)
(*Export the result as animated GIF.*)
Print@Export[
FileNameJoin[{"mackay_gaussian_2006-bi-prior.gif"}],
frames,
AnimationRepetitions->Infinity,
"DisplayDurations"->.1
]
Raw
mackay_gaussian_2006-multi-posterior.wls
#!/usr/bin/env wolframscript
(* ::Package:: *)
(* ::Text:: *)
(*This script generates the fourth animated figure on my blog post "mackay_gaussian_2006".*)
Print["https://jem-mosig.com/2019/03/mackay_gaussian_2006/"]
(* ::Text:: *)
(*Generate all frames of the animation, one sample from the bi-variate normal distribution for each frame.*)
Print["Generating figure. This may take a minute..."]
frames = BlockRandom[
SeedRandom[2019];
With[{x=Delete[Subdivide[-5,5,40],13],\[Sigma]=1,l=1,x0=-2,y0=-1},
Module[{\[Mu],K,Kinv,A1,A2,B,\[Mu]post,Kpost},
\[Mu]=ConstantArray[0,Length[x]];K=N[Outer[\[Sigma]^2 Exp[-(#1-#2)^2/(2l^2)]&,Join[{x0},x],Join[{x0},x]],40];
Kinv=Inverse[K];
A1=Kinv[[1;;1,1;;1]];
B=Kinv[[1;;1,2;;]];
A2=Kinv[[2;;,2;;]];
Assert[
AllTrue[Flatten@Chop[Kinv-ArrayFlatten[{{A1,B},{Transpose[B],A2}}]],#==0&]
];
\[Mu]post=-Inverse[A2].Transpose[B].{{y0}}//Flatten;
Kpost=Inverse[A2];
Quiet[
With[{dist=MultinormalDistribution[\[Mu]post,Kpost],stddiv=Sqrt@Diagonal[Kpost]},
Table[
With[{y=RandomVariate[dist]},
Show[{
ListPlot[
Transpose[{x,y}],
Joined->True,
InterpolationOrder->2,
PlotStyle->Thick,
PlotRange->{-3,3},
Axes->False,Frame->True,
Mesh->Full,
MeshStyle->Directive[ColorData[112,2],PointSize[Large]],
Epilog->{PointSize[Large],Point[{x0,y0}]},
FrameLabel->{"x","f \[Distributed] P(y|\!\(\*SubscriptBox[\(y\), \(0\)]\),\[CapitalSigma])"}
],
ListPlot[{
SortBy[Append[Transpose[{x,\[Mu]post+stddiv}],{x0,y0}],First],
SortBy[Append[Transpose[{x,\[Mu]post-stddiv}],{x0,y0}],First],
Transpose[{x,\[Mu]post}]
},
Joined->True,
PlotStyle->{Lighter@ColorData[112,2],Lighter@ColorData[112,2],Directive[Dashed,ColorData[112,2]]},
Filling->{1->{2}}
]
}]
],
100
]
],
{MultinormalDistribution::posdefprm,Transpose::nmtx}
]
]
]
]
(* ::Text:: *)
(*Sort the frames to appear visually more pleasing.*)
ordering = Quiet[
Last@FindShortestTour@Table[
FirstCase[frame,Line[pts_]:>Interpolation[pts][Subdivide[-4.5,4.5,10]],None,Infinity],
{frame,frames}
],
{InterpolatingFunction::dmval}
];
frames = frames[[ordering]];
(* ::Text:: *)
(*Export the result as animated GIF.*)
Print@Export[
FileNameJoin[{"mackay_gaussian_2006-multi-posterior.gif"}],
frames,
AnimationRepetitions->Infinity,
"DisplayDurations"->.5
]
Raw
mackay_gaussian_2006-multi-prior.wls
#!/usr/bin/env wolframscript
(* ::Package:: *)
(* ::Text:: *)
(*This script generates the third animated figure on my blog post "mackay_gaussian_2006".*)
Print["https://jem-mosig.com/2019/03/mackay_gaussian_2006/"]
(* ::Text:: *)
(*Generate all frames of the animation, one sample from the bi-variate normal distribution for each frame.*)
Print["Generating figure. This may take a minute..."]
frames = BlockRandom[
SeedRandom[2019];
With[{x=Subdivide[-5,5,40],\[Sigma]=1,l=1},
Module[{\[Mu],K,Kinv,A1,A2,B,\[Mu]post,Kpost},
\[Mu]=ConstantArray[0,Length[x]];K=N[Outer[\[Sigma]^2 Exp[-(#1-#2)^2/(2l^2)]&,x,x],40];
Quiet[
With[{dist=MultinormalDistribution[\[Mu],K],stddiv=Sqrt@Diagonal[K]},
Table[
With[{y=RandomVariate[dist]},
Show[{
ListPlot[
Transpose[{x,y}],
Joined->True,
InterpolationOrder->1,
PlotStyle->Thick,
PlotRange->{-3,3},
Axes->False,Frame->True,
Epilog->{
ColorData[112,2],
PointSize[Large],
Point[Transpose[{x,y}]]
},
FrameLabel->{"x","f \[Distributed] P(y|\[CapitalSigma])"}
],
ListPlot[{
Transpose[{x,\[Mu]+stddiv}],
Transpose[{x,\[Mu]-stddiv}],
Transpose[{x,\[Mu]}]
},
Joined->True,
PlotStyle->{Lighter@ColorData[112,2],Lighter@ColorData[112,2],Directive[Dashed,ColorData[112,2]]},
Filling->{1->{2}}
]
}]
],
100
]
],
{MultinormalDistribution::posdefprm,Transpose::nmtx}
]
]
]
];
(* ::Text:: *)
(*Sort the frames to appear visually more pleasing.*)
ordering = Quiet[
Last@FindShortestTour@Table[
FirstCase[frame,Line[pts_]:>Interpolation[pts][Subdivide[-4.5,4.5,10]],None,Infinity],
{frame,frames}
]
];
frames=frames[[ordering]];
(* ::Text:: *)
(*Export the result as animated GIF.*)
Print@Export[
FileNameJoin[{"mackay_gaussian_2006-multi-prior.gif"}],
frames,
AnimationRepetitions->Infinity,
"DisplayDurations"->.5
]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment