<<MaTeX`

MaTeX["CDF(x):=\\int_0^x pdf(z)dz=\\xi",FontSize->24]
MaTeX["CDF(x)=\\xi\\quad  \\Rightarrow \\quad \\int_0^x pdf(z)dz=\\xi ",FontSize->24]

The pdf can be seen as the Jacobian factor from the original variable z and the uniformly distributed variable ξ.

Define a PDF (with known integral)

pdf[λ_,z_]=λ*Exp[-λ*z]

Define two probability density function from the same non-normalized “pdf”, one with its normalization factor given by the analytic integral, computed with Integrate[], the other with a numerical integrator NIntegrate[].

Clear[{npdf,Nnpdf}]
npdf[λ_,x_]=pdf[λ, x]/Integrate[pdf[λ, z], {z, 0, 1}]
Nnpdf[λ_,x_]:=pdf[λ, x]/NIntegrate[pdf[λ, z], {z, 0, 1}]

In this case there is a useful difference between the `SetDelayed` operation `:=` and the ordinary `Set` operation `=`
(HyperlinkURL["http://reference.wolfram.com/language/ref/SetDelayed.html", http://reference.wolfram.com/language/ref/SetDelayed.html])
(HyperlinkURL["http://reference.wolfram.com/language/ref/Set.html", http://reference.wolfram.com/language/ref/Set.html])

Define two Cumulative functions, one with its normalization factor given by the analytic integral, computed with Integrate[], the other with a numerical integrator NIntegrate[].

Clear[cdf]
cdf[λ_,x_]=Integrate[npdf[λ,z],{z,0,x}]
Clear[Ncdf]
Ncdf[λ_,x_Real]:=NIntegrate[npdf[λ,z],{z,0,x}]

npdf[1,1.1]
Nnpdf[1,1.1]

cdf[1,1.1]
Ncdf[1,1.1]

cdf[z,1]

cdf[z,1.1]

cdf[z,q]

Display the CDF for different values of λ

lambdas={1,2,3,10}
Plot[Evaluate[Table[Ncdf[λ,w],{λ,lambdas}]],{w,0,1},PlotRange->All,PlotLegends->lambdas,Frame->True,FrameLabel->{"x","CDF(λ,x)"}]
Unset[lambdas]

uniformValues=RandomReal[{0,1},3000];

Use FindRoot as method to solve numerically the equation for ξ=CDF(x)

?FindRoot

FindRoot[Ncdf[1,x]=={0.2,0.22},{x,0.5}]

?RealValuedNumberQ

The solution might be in principle complex. Indeed using NSolve instead of FindRoot it would return several unwanted complex solution. We make a filter to remove them, just in case.

Select[{x/.FindRoot[cdf[1,x]==0.2,{x,0.5}]},RealValuedNumberQ]

Select[{x/.FindRoot[Ncdf[1,x]==0.2,{x,0.5}]},RealValuedNumberQ]

Select[x/.NSolve[cdf[1,x]==0.2,x],RealValuedNumberQ]

?AbsoluteTiming

We can time the execution using the numerically normalized CDF and the analytically normalized CDF. The latter is 80 times faster.

AbsoluteTiming[psols=ParallelMap[(Select[{x/.FindRoot[cdf[1,x]==#,{x,0.5}]},RealValuedNumberQ])&,uniformValues]]

AbsoluteTiming[psols=ParallelMap[(Select[{x/.FindRoot[Ncdf[1,x]==#,{x,0.5}]},RealValuedNumberQ])&,uniformValues]]

Dimensions[psols]
Dimensions[Flatten[psols]]

The values obtained with the Jacobian method follow the distribution we wanted!

observations= Flatten[psols];
Show[Histogram[observations,{0.1},"PDF"],
Plot[(*Length[observations]**)npdf[1,x],{x,0,1}]]

observations= Flatten[psols];
Show[Histogram[observations,{0.1},"PDF"],
Plot[(*Length[observations]**)npdf[1,x],{x,0,1}]]

Oscillating PDF: Accept-Reject method

When the integral of the PDF is not easy to obtain or we want to not use numerical solutions like FindRoot or NSolve we can use the Accept-Reject method.

Find a constant C that is larger than any value taken by  the pdf in its domain, which we take to be [0,1]. 
Generate (Subscript[N, p]) values (Subscript[x, i]) in the domain, which we take to be [0,1].
Generate (Subscript[N, p]) uniform values (Subscript[u, i]) in the range [0,1]
Accept only the x values in the set {(Subscript[x, i])} for which pdf((Subscript[x, i])) is greater than C (Subscript[u, i]), scanning all i in [1,(Subscript[N, p])].

MaTeX["pdf(x_i)>C \\cdot u_i",FontSize->33]

OscillatingPdf[λ_,z_]=1+Cos[-λ*z]
Integrate[OscillatingPdf[λ,z],{z,0,1}]
NormalizedOscillatingPdf[λ_,z_]=(1+Cos[-λ*z])/Integrate[OscillatingPdf[λ,z],{z,0,1}]

xValues=RandomReal[{0,1},10^5];
ξValues=RandomReal[{0,1},10^5];

pdfValues=NormalizedOscillatingPdf[4,#]&@xValues
Const=Max[pdfValues]

test=Map[#>0&,pdfValues-Const*ξValues];
filtered=Select[Transpose[{xValues,test}],(#[[2]]==True)&];
Length[filtered]/Length[pdfValues]//N
Show[Histogram[filtered[[All,1]],{0.0025},"PDF"],Plot[NormalizedOscillatingPdf[4,z],{z,0,3}]]

Stratified sampling to populate the tails more accurately

xValues1=Select[xValues,(#<=0.5)&]
xValues2=Select[xValues,(#>0.5)&]
pdfValues1=NormalizedOscillatingPdf[4,#]&@xValues1
pdfValues2=NormalizedOscillatingPdf[4,#]&@xValues2
Const1=Max[pdfValues1]
Const2=Max[pdfValues2]

test1=Map[#>0&,pdfValues1-Const1*ξValues[[1;;Length[pdfValues1]]]];
test2=Map[#>0&,pdfValues2-Const2*ξValues[[Length[pdfValues1]+1;;-1]]];
filtered1=Select[Transpose[{xValues1,test1}],(#[[2]]==True)&];
filtered2=Select[Transpose[{xValues2,test2}],(#[[2]]==True)&];
Length[filtered1]/Length[pdfValues1]//N
Length[filtered2]/Length[pdfValues2]//N
totalCounts=Length[filtered1[[All,1]]]+Length[filtered2[[All,2]]]
Show[Histogram[filtered1[[All,1]],{0.025},"Count"],Histogram[filtered2[[All,1]],{0.025},"Count"],(*Plot[(totalCounts)*NormalizedOscillatingPdf[4,z],{z,0,1}],*)PlotRange->{{0,1},All}]

The heights of the two sub-regions do not make the histogram continuous!
We need to weigh the data from the two sub-regions so that the histogram is continuous and it preserves the area of the function.  The function `WeightedData[]` is our friend.

data={8,3,5,4,9,0,4,2,2,3};
w={0.15,0.09,0.12,0.10,0.16,0.,0.11,0.08,0.08,0.09};
wd=WeightedData[data,w]
Mean[data]
Mean[wd]

test1=Map[#>0&,pdfValues1-Const1*ξValues[[1;;Length[pdfValues1]]]];
test2=Map[#>0&,pdfValues2-Const2*ξValues[[Length[pdfValues1]+1;;-1]]];
filtered1=Select[Transpose[{xValues1,test1}],(#[[2]]==True)&];
filtered2=Select[Transpose[{xValues2,test2}],(#[[2]]==True)&];
Length[filtered1]/Length[pdfValues1]//N
Length[filtered2]/Length[pdfValues2]//N
totalCounts=Length[filtered1[[All,1]]]+Length[filtered2[[All,2]]]
overall=22
wd1=WeightedData[filtered1[[All,1]],Table[overall*Const1/(Const2+Const1),Length[filtered1[[All,1]]]]];
wd2=WeightedData[filtered2[[All,1]],Table[overall*Const2/(Const2+Const1),Length[filtered2[[All,1]]]]];
binsize=0.0025;
Show[
Histogram[wd1,{binsize},"Count"],Histogram[wd2,{binsize},"Count",ChartStyle->Hue[.5]],(*Histogram[filtered[[All,1]],{0.025},"Count",ChartStyle->Opacity[0.1,Hue[.5]]],*)Plot[(Total[wd1["Weights"]]+Total[wd2["Weights"]])/(1/binsize)*NormalizedOscillatingPdf[4,z],{z,0,1},PlotStyle->Dotted],PlotRange->{{0,1},All}]

binsize=1/400.;test1=Map[#>0&,pdfValues1-Const1*ξValues[[1;;Length[pdfValues1]]]];
test2=Map[#>0&,pdfValues2-Const2*ξValues[[Length[pdfValues1]+1;;-1]]];
filtered1=Select[Transpose[{xValues1,test1}],(#[[2]]==True)&];
filtered2=Select[Transpose[{xValues2,test2}],(#[[2]]==True)&];
Length[filtered1]/Length[pdfValues1]//N
Length[filtered2]/Length[pdfValues2]//N
totalCounts=Length[filtered1[[All,1]]]+Length[filtered2[[All,2]]]
l=Length[filtered[[All,1]]];
l1=Length[filtered1[[All,1]]];
l2=Length[filtered2[[All,1]]];
wd1=WeightedData[filtered1[[All,1]],Table[Const1/(Const2+Const1)*(l1+l2),l1]];
wd2=WeightedData[filtered2[[All,1]],Table[Const2/(Const2+Const1)*(l1+l2),l2]];
wd=WeightedData[filtered[[All,1]],Table[l,l]];
Show[
Histogram[wd1,{binsize},"Count"],
Histogram[wd2,{binsize},"Count",ChartStyle->Hue[.5]],Histogram[wd,{binsize},"Count",ChartStyle->Opacity[0.1,Hue[.5]]],Plot[(Total[wd1["Weights"]]+Total[wd2["Weights"]])/(1/binsize)*NormalizedOscillatingPdf[4,z],{z,0,1}],
Plot[(Total[wd["Weights"]])/(1/binsize)*NormalizedOscillatingPdf[4,z],{z,0,1}],PlotRange->{{0,1},All}]