Algorithm conversion to Matlab Code

20 次查看(过去 30 天)
Mar One
Mar One 2022-5-21
回答: Todd 2025-9-14
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

回答(4 个)

Steven Lord
Steven Lord 2022-5-21
See this post on Cleve Moler's blog.
  2 个评论
Mar One
Mar One 2022-5-21
i need the matlab code for the algorithm if posible
Walter Roberson
Walter Roberson 2022-5-22
Calling randn() is the matlab code that implements the algorithm.

请先登录,再进行评论。


Walter Roberson
Walter Roberson 2022-5-22

Soufi
Soufi 2024-11-9
how can i traduct algorithme to matlab

Todd
Todd 2025-9-14
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={17}
}
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

类别

Help CenterFile Exchange 中查找有关 Random Number Generation 的更多信息

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by