vervaatsim

vervaatsim returns a Vervaat perpetuity.

expand all in page

Syntax

Description

This function allows to simulate exactly from a Vervaat perpetuity distribution with parameter $\beta \in (0,+\infty)$. The simulation method, by K. Cloud, M. Huber (2018), is accurate and very efficient: it is $\mathcal{O}(\beta \log(\beta))$, in the sense that it uses $T$ uniform random variates, where $\mathbf{E}[T] \le \mathcal{O}(\beta \log(\beta))$. Therefore the algorithm is good even for $\beta \gg 1$.

Remark: this function is recursive and cannot be easily extended to the generation of several vervaat numbers.

The FSDA code has been translated from the R version of the authors.

Below we give briefly some definitions, background and motivations. In the references, we also include several key works on the distribution.

example

y =vervaatsim(betav) Density of Vervaat perpetuity with $\beta=1$ and no optional parameters.

example

y =vervaatsim(betav, steps) Density of Vervaat perpetuity with $\beta=1$ and 10 chain steps.

example

y =vervaatsim(betav, steps, d) Density of Vervaat perpetuity with $\beta=1$, 10 chain steps and initial chain value d=2.

Examples

expand all

  • Density of Vervaat perpetuity with $\beta=1$ and no optional parameters.

    • clear all;
close all;
betav = 1;
y=vervaatsim(betav)

  • Density of Vervaat perpetuity with $\beta=1$ and 10 chain steps.

    • clear all;
close all;
betav = 1;
steps = 10;
y=vervaatsim(betav,steps)

  • Density of Vervaat perpetuity with $\beta=1$, 10 chain steps and initial chain value d=2.

    • clear all;
close all;
betav = 1;
steps = 10;
d = 2;
y=vervaatsim(betav,steps,d)

    Related Examples

    expand all

  • N=5000 random values extracted from two Vervaat perpetuities.

    • Vervaat parameters are: $\beta = 1$ and $\beta = 10$. 

    % The superimposed normal kernel density is just for illustration.
% The same for the superimposed cdf, made using FSDA function exactcdf.
clear all; 
close all;
betav10 = 10;
betav01 = 1;
N = 5000;
y10 = zeros(N,1); y01 = y10;
for i=1:N
y10(i) = vervaatsim(betav10);
y01(i) = vervaatsim(betav01);
end
pdy10 = fitdist(y10(:),'Kernel','Kernel','normal','Support','positive');
x10   = (1:N)*(5*betav10)/N;
pdf10 = pdf(pdy10,x10);
pdy01 = fitdist(y01(:),'Kernel','Kernel','normal','Support','positive');
x01   = (1:N)*(5*betav01)/N;
pdf01 = pdf(pdy01,x01);
figure;
h1=subplot(2,1,1);
plot(x10,pdf10,'-','LineWidth',2); ylim([0,1]);
% now superimpose the cdf using function exactcdf
cdf10 = exactcdf(x10,y10);
hold on;
plot(x10,cdf10,'-r','LineWidth',2); ylim([0,1]);
hold off;
% superimpose the histogram
if ~verLessThan('matlab','1.7.0')
hold on
h=histogram(y10);
h.Normalization='pdf';
h.BinWidth=0.02;
h.EdgeColor='none'; 
hold off
end
drawnow;
h2=subplot(2,1,2);
plot(x01,pdf01,'-','LineWidth',2); ylim([0,1]);
% now superimpose the cdf using function exactcdf
cdf01 = exactcdf(x01,y01);
hold on;
plot(x01,cdf01,'-r','LineWidth',2); ylim([0,1]);
hold off;
% superimpose the histogram
if ~verLessThan('matlab','1.7.0')
hold on
h=histogram(y01);
h.Normalization='pdf';
h.BinWidth=0.02;
h.EdgeColor='none'; 
hold off
end
drawnow;
title(h1,'$$\beta = 10$$ with cdf superimposed' ,'Fontsize',20,'interpreter','latex');
title(h2,'$$\beta = 1 $$ (Dickman) with cdf superimposed'  ,'Fontsize',20,'interpreter','latex');
    Click here for the graphical output of this example (link to Ro.S.A. website)

  • A rough time test on the method.

    • clear all;
close all;
betav = 1;
n = 100;
nmax    = 500;
betamax = 10;
rng('default')
rng(123);
N     = randi(nmax,n,1);
betav = randi(betamax,n,1);
ti = tic;
for i=1:n;
for j=1:N(i)
y1(j) = vervaatsim(betav(i));
end
end
tf = toc(ti);
disp(['Cloud-Huber: etime = ' num2str(tf)]);

  • Density of Vervaat perpetuity: comparison with original R function.

    • clear all;
close all;
% set common parameters
betav = 1;
steps = 10;
d = 3;
% ensure same random numbers in MATLAB and R 
rn = mtR(1);    
% Now compute the veervaat in MATLAB
y=vervaatsim(betav,steps,d)
% And finally compute the vervaat in R
disp('To verify coherence of results with the R implementation,');
disp('execute in R the following four lines,'); 
disp('which allow to replicate the same random numbers;');
disp('then execute the vervaat after sourcing the following R code: ') ;
disp(' ');
disp('NGkind("Mersenne-Twister") #set "Mersenne-Twister" "Inversion"');
disp('set.seed(0) ');
disp('state = .Random.seed ');
disp('runif(1) ');
disp('vervaat(beta = 1,steps = 10,d = 3) ');
disp(' ');
    y =

    1.8183

To verify coherence of results with the R implementation,
execute in R the following four lines,
which allow to replicate the same random numbers;
then execute the vervaat after sourcing the following R code: 
 
NGkind("Mersenne-Twister") #set "Mersenne-Twister" "Inversion"
set.seed(0) 
state = .Random.seed 
runif(1) 
vervaat(beta = 1,steps = 10,d = 3)

    Input Arguments

    expand all

    betav — Distribution parameter value. Positive integer.

    The parameter of the Vervaat family. Dickman function is obtained for betav = 1.

    Data Types: single| double

    Optional Arguments

    steps — Markov chain step. Positive integer.

    The inital number of steps to run to move the chain forward. Default is steps = 1.

    Example: steps=5.

    Data Types: Scalar.

    d — Chain value at time 0. Positive integer or -1.

    The value of the dominating chain at time 0. Default is $D_0 \leftarrow x_0−1+G$, where $G \sim \mbox{Geo}(1/2)$ and $x_0 = \frac{1 + (2/3)^{1/\beta}}{1 - (2/3)^{1/\beta}}$ ($\mbox{Geo}$ represents the geometric distribution).

    Example: d=10.

    Data Types: Scalar.

    Output Arguments

    expand all

    y —The Vervaat perpetuity value. Scalar

    The estimated distribution value. Data Types - Double.

    More About

    expand all

    Additional Details

    A perpetuity is a random variable of the form: \[ Y = W_1 + W_1 \cdot W_2 + W_1 \cdot W_2 \cdot W_3 + \ldots \]

    where the $W_i$ are an independent, identically distributed (iid) sequence of random variables. If each $W_i$ has the same distribution, say $W_i \sim W$, then $Y \sim W(1 + Y)$ for $Y$ and $W$ independent.

    We are interested in this distribution because the running time of the Quickselect algorithm of Hoare, for finding the order statistics in a numerical array, approaches asymptotically a particular perpetuity, called Dickman distribution, with $W \sim Unif([0,1])$. Unfortunately such distribution has no closed form.

    The Dickman distribution can be also seen as a special case of Vervaat perpetuity, which is such that $W_i \sim U^{1/\beta}$ for some $\beta \in (0,\infty)$ for $U \sim Unif([0,1])$. In other words, the Dickman distributon is a Vervaat perpetutiy with $\beta = 1$.

    A generalization of the perpetuity takes the form

    \[ Y = \sum_{n=0}^{\infty}A_n \prod_{i}^{n} W_i\]

    with $A_n$ not necessarily equal to $1$, which is known as Takacs distribution.

    References

    Cloud, K. and Huber, M. (2018), Fast Perfect Simulation of Vervaat Perpetuities, "Journal of Complexity", Vol. 42, pp. 19-30.

    Devroye, L. (2001), Simulating perpetuities, "Methodology And Computing In Applied Probability", Vol. 3, Num. 1, pp. 97–115.

    Fill, J. A. and Huber, M. (2010), Perfect simulation of Vervaat perpetuities, "Electronic Journal of Probability", Vol. 15, pp. 96–109.

    Devroye, L. and Fawzi, O. (2010), Simulating the Dickman distribution, "Statistics and Probability Letters", Vol. 80, pp. 242–247.

    Blanchet, J. H. and Sigman, K. (2011), On exact sampling of stochastic perpetuities, "Journal of Applied Probability", Vol. 48A, pp. 165–182.

    Takacs, L. (1955), On stochastic processes connected with certain physical recording apparatuses. "Acta Mathematica Academiae Scientificarum Hungarica", Vol. 6, pp 363-379.

    Barabesi, L. and Pratelli, L. (2019), On the properties of a Takacs distribution, "Statistics and Probability Letters", Vol. 148, pp. 66-73.

    See Also

    | |

    This page has been automatically generated by our routine publishFS