Instantly share code, notes, and snippets.

# JEM-Mosig/blog-mackay_gaussian_2006.md

Last active April 2, 2019 10:26
Show Gist options
• Save JEM-Mosig/39ac1a9f88d9bd7bc18c168aeb30dca1 to your computer and use it in GitHub Desktop.
blog-mackay_gaussian_2006

# 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.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #!/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 ]
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #!/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 ]
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #!/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 ]
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #!/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 ]