Algorithm conversion to Matlab Code

Hello, i have this ziggurat algorithm that i wish to convert to a matlab code if someone would know the correct matlab code for it

Answer 1

[Steven Lord]

See this post on Cleve Moler's blog.

[Mar One]

i need the matlab code for the algorithm if posible

[Walter Roberson]

Calling randn() is the matlab code that implements the algorithm.

Answer 2

[Walter Roberson]

https://people.sc.fsu.edu/~jburkardt/m_src/ziggurat/ziggurat.html

Answer 3

[Soufi]

how can i traduct algorithme to matlab

[Walter Roberson]

Source code is at https://people.sc.fsu.edu/~jburkardt/m_src/ziggurat/ziggurat.html

Answer 4

[Todd]

Based on wikipedia, this is most likely from

@article{JSSv005i08,
 title={The Ziggurat Method for Generating Random Variables},
 volume={5},
 url={https://www.jstatsoft.org/index.php/jss/article/view/v005i08},
 doi={10.18637/jss.v005.i08},
 abstract={We provide a new version of our ziggurat method for generating a random variable from a given decreasing density. It is faster and simpler than the original, and will produce, for example, normal or exponential variates at the rate of 15 million per second with a C version on a 400MHz PC. It uses two tables, integers k<sub>i</sub>, and reals w<sub>i</sub>. Some 99% of the time, the required x is produced by: Generate a random 32-bit integer j and let i be the index formed from the rightmost 8 bits of j. If j < k, return x = j x w<sub>i</sub>. We illustrate with C code that provides for inline generation of both normal and exponential variables, with a short procedure for settting up the necessary tables.},
 number={8},
 journal={Journal of Statistical Software},
 author={Marsaglia, George and Tsang, Wai Wan},
 year={2000},
 pages={1–7}
}

The algorithm's source code is then

\usepackage{algorithm}
\usepackage{algpseudocode}

\begin{algorithm}
\caption{Original Ziggurat Algorithm (optimizations omitted for clarity)}
\RequireItem{$N$}{$\triangleright$ Number of regions}
\RequireItem{$x[0..N],\,y[0..N]$}{$\triangleright$ Coordinates of regions}
\RequireItem{$\mathrm{PDF}(\cdot)$}{$\triangleright$ Probability density function}
\RequireItem{$\mathrm{RandReal}()$}{$\triangleright$ A uniform real in $[0,1]$}
\RequireItem{$\mathrm{RandInteger}(N)$}{$\triangleright$ A uniform integer in $[0,N]$}

\begin{algorithmic}[1]
\Function{Ziggurat}{}
  \Loop
    \State $j \gets \mathrm{RandInteger}(N)$
    \State $x \gets x[j] \times \mathrm{RandReal}()$
    \If{$x < x[j+1]$} 
      \State \Return $x$
    \ElsIf{$j \neq 0\;$ and $\;\mathrm{RandReal}() \times (y[j+1] - y[j]) < \mathrm{PDF}(x) - y[j]$}
      \State \Return $x$
    \ElsIf{$j = 0$}
      \State \Return \Call{Tail}{$x[1]$}
    \EndIf
  \EndLoop
\EndFunction
\end{algorithmic}
\end{algorithm}

with Matlab code:

function z = zigguratSample(N, x, y, pdf_func)
% ZIGGURATSAMPLE   Sample one variate using the Ziggurat algorithm
%   Inputs:
%     N        : number of regions (integer)
%     x(0:N)   : right-edge coordinate of regions
%     y(0:N)   : heights of regions
%     pdf_func : function handle, pdf_func(x) returns f(x)
%
%   Output:
%     z        : a variate drawn from the target density

    while true
        j = randi([0, N]);          % uniform integer in 0..N
        u = rand();                 % uniform real in [0,1)
        xj = x(j+1);                % MATLAB indices are 1-based; x(1) = x[0], etc
        xjp1 = x(j+2);              % x[j+1]
        yj = y(j+1);
        yjp1 = y(j+2);

        x_star = xj * rand();       % propose horizontal coordinate

        if x_star < xjp1
            z = x_star;
            return;
        elseif j ~= 0 && rand() * (yjp1 - yj) < (pdf_func(x_star) - yj)
            z = x_star;
            return;
        elseif j == 0
            z = tailSample(x(2), pdf_func);  % tail case, using x[1] in 0-based
            return;
        end
        % else continue loop
    end
end


function f = pdf_normal(x)
    % Example: standard normal pdf
    f = exp(-0.5 * x.^2) / sqrt(2*pi);
end

function z = tailSample(x1, pdf_func)
    % Tail sampler; for normal, use exponentials etc.
    while true
        u = rand();
        v = rand();
        x = -log(u) / x1;
        y = -log(v);
        if 2*y >= x^2
            z = x + x1;
            return;
        end
    end
end