Code covered by the BSD License

# Random numbers from a user defined piecewise polynomial distribution

### Luigi Sanguigno (view profile)

RANDPDFS returns pseudorandom values drawn from an arbitrary distribution defined by the pp-form.

[r,varargout]=randpdfs(ppdf,varargin)
```function [r,varargout]=randpdfs(ppdf,varargin)
% RANDPDFS Pseudorandom numbers from an arbitrary distribution.
%     R = RANDPDFS(PPDF,N) returns an N-by-N matrix containing pseudorandom values drawn
%     from an arbitrary distribution defined by the pp-form PPDF. RANDPDFS(PPDFM,N)
%     or RANDPDFS(PPDF,[M,N]) returns an M-by-N matrix.  RAND(M,N,P,...) or
%     RANDPDFS(PPDF,[M,N,P,...]) returns an M-by-N-by-P-by-... array.
%     RANDPDFS(PPDF,SIZE(A)) returns an array the same size as A.
%     [R,RANDPP] = RANDPDFS(PPDF,...) provides the pp-form used to pass
%     from the uniform distribution to the user defined one.
%
%     Note: The size inputs M, N, P, ... should be nonnegative integers.
%     Negative integers are treated as 0.
%
%     R = RANDPDFS(..., 'double') or R = RANDPDFS(..., 'single') returns an array of
%     uniform values of the specified class.
%
%     Compatibility Note: RANDPDFS uses FNINT to compute the cumulative
%     distribution function, therefore spline toolbox has to be accessible
%     to allow RANDPDFS working properly.
%

% Compute cumulative distribution
ppcdf=fnint(ppdf);
ppcdf.coefs=ppcdf.coefs/ppval(ppcdf,ppcdf.breaks(end));

% Calculate the cumulative distribution function inverse
invppcdf=spline(ppval(ppcdf,ppcdf.breaks),...
cat(2,1/ppval(ppdf,ppcdf.breaks(1)),ppcdf.breaks,1/ppval(ppdf,ppcdf.breaks(end))));

% Compute random sequence
r=ppval(invppcdf,rand(varargin{:}));

if nargout==2
varargout{1}=invppcdf;
end```