# How to Sample from a Probability Density Function

I've done a pretty good job working my way through and understanding the various ways to calculate the differences between probability distributions through Kullback Leibler Divergence and was happy to get the same results with my example in Mathematica 11.3.

distP = NormalDistribution[];
distQ = GumbelDistribution[];
p[x_]:=PDF[distP, x]
q[x_]:=PDF[distQ, x]

(* this is the difference between distributions *)
solA = Expectation[Log[p[x]/q[x]],x\[Distributed]distP]
(* Same as above *)
solB = Integrate[p[x] Log[p[x]/q[x]],{x,-Infinity,Infinity}]
solA == solB
(* True *) However, when I came upon the following image near the bottom of the page and tried to replicate it in my notebook, I found I was a bit lost as to how to sample points from a probability distribution in the form of p(x)? How to sum and average the points makes sense, but not how to sample them N times sequentially as required by the problem to converge on the difference. • One of the problems with the blog you reference is that the statement with the limit is nonsensical. The left-hand side is about an average of a specific set of realizations (or if being charitable a statement about an infinite number of random variables or if not being charitable "not well defined") and the right-hand side is a well-defined expectation a single random variable.Totally different critters and not at all equivalent. The graph gives the impression that the expectation is off unless you have a large number of samples:.not true. Learn Statistics from a Statistician not an engineer. – JimB Apr 15 at 3:02
• Does TransformedDistribution[Log[p[x]/q[x]], x \[Distributed] distP] work for you? You should be able to use RandomVariate on that. – Sjoerd Smit Apr 15 at 13:21
• @Sjoerd Smit I had never seen that function before but it indeed starts to approach as n increases. RandomVariate[TransformedDistribution[Log10[p[x]/q[x]],x\[Distributed]distP],10^5] // Mean – BBirdsell Apr 16 at 4:54

Here is one way to do it.

SeedRandom@2;
distP = NormalDistribution[];
distQ = GumbelDistribution[];
p[x_] := PDF[distP, x]
q[x_] := PDF[distQ, x]
x = Table[RandomVariate[distP, n], {n, 1, 10000, 100}];
logLR = Mean /@ Log10[p[x]/q[x]];
nValues = Log10@Range[1, 10000, 100];
ListLinePlot[Transpose@{nValues, logLR}, PlotRange -> {0, 0.5}, Frame -> True,
PlotTheme -> "Detailed"] Or use

ListLogLinearPlot[logLR, PlotRange -> {0, 0.5}, Joined -> True,
Frame -> True, PlotTheme -> "Detailed"]


Here is the confirmation of website's example.

SeedRandom@2;
p[x_] := PDF[distP, x]
q[x_] := PDF[distQ, x]
x = Table[RandomVariate[distP, n], {n, 1, 10000, 100}];
logLR = Mean /@ Log10[p[x]/q[x]];
nValues = Log10@Range[1, 10000, 100];
ListLinePlot[Transpose@{nValues, logLR}, PlotRange -> {0, 0.5}, Frame -> True,
PlotTheme -> "Detailed"] • I didn't know that a list structured like x's could just be put into those p/q functions and they would just map over them predictably. the Log10 also appears to be important. Thanks for taking the time to write this up. – BBirdsell Apr 15 at 4:17
• See my edit. It is more consistent with websites example.. – OkkesDulgerci Apr 15 at 12:52

Not sure why I'm getting something slightly different from the previous answer. Is it supposed to converge at solA, or did I misunderstand the question?

distP = NormalDistribution[0, 1];
distQ = GumbelDistribution[0, 1];
p[x_] := PDF[distP, x]
q[x_] := PDF[distQ, x]
solA = Expectation[Log[p[x]/q[x]], x \[Distributed] distP];

logLR[s_] := Module[
{sample = RandomVariate[distP, s]},
Sum[Log[p[sample[[i]]]/q[sample[[i]]]], {i, s}]/s
]

l = Table[{Log[10, s], logLR[s]}, {s, 1, 3000, 10}];

Show[
ListPlot[l, Joined -> True, PlotRange -> {0, .7}],
Graphics[{Thick, Dashed, Line[{{0, N[solA]}, {3.5, N[solB]}}]}]] • It's suppose to converge to the Difference as far as I know. Needed Log10 for everything to work out for me in the end. – BBirdsell Apr 15 at 4:43