randreject_yp

 Calculate N times the random variable x whose p.d.f. is f(x). f(x) may
 not be normalized to unity on the integration domain [xmin,xmax], where
 xmin < xmax.

 Input:

   - fun_handle : handle to the function f(x) which may be parametrized and must be expressed as f(a1,a2,a3,...,[xmin,xmax],x). It is better that f(x) is normalized to unity on the integration domain [xmin,xmax].
   - args_in : N-1 first arguments of the function f(a1,a2,a3,...,[xmin,xmax],x): a1,a2,a3,...,[xmin,xmax]. The last arguments of args_in must be the range of the random variable [xmin,xmax] on which f(x) is calculated.
   - N : number of trials [1,1]
     (default = 100000)
   - display_mode : display the numerical and theoretical p.d.f and display the distribution is > 0 and is the number of bin of the histogram (0: no, > 1, yes) [1,1]
     (default = 0)

 Output:

   - y : random variable [N,1]
   - y_th : theoretical random variable [N,1]
   - py_th : theoretical p.d.f. normalized to 1 [N,1]

   by Y. Peysson (CEA/DSM/IRFM) (yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM) (joan.decker@cea.fr)

0001 function [y,y_th,py_th] = randreject_yp(fun_handle,args_in,N,display_mode)
0002 %
0003 % Calculate N times the random variable x whose p.d.f. is f(x). f(x) may
0004 % not be normalized to unity on the integration domain [xmin,xmax], where
0005 % xmin < xmax.
0006 %
0007 % Input:
0008 %
0009 %   - fun_handle : handle to the function f(x) which may be parametrized and must be expressed as f(a1,a2,a3,...,[xmin,xmax],x). It is better that f(x) is normalized to unity on the integration domain [xmin,xmax].
0010 %   - args_in : N-1 first arguments of the function f(a1,a2,a3,...,[xmin,xmax],x): a1,a2,a3,...,[xmin,xmax]. The last arguments of args_in must be the range of the random variable [xmin,xmax] on which f(x) is calculated.
0011 %   - N : number of trials [1,1]
0012 %     (default = 100000)
0013 %   - display_mode : display the numerical and theoretical p.d.f and display the distribution is > 0 and is the number of bin of the histogram (0: no, > 1, yes) [1,1]
0014 %     (default = 0)
0015 %
0016 % Output:
0017 %
0018 %   - y : random variable [N,1]
0019 %   - y_th : theoretical random variable [N,1]
0020 %   - py_th : theoretical p.d.f. normalized to 1 [N,1]
0021 %
0022 %   by Y. Peysson (CEA/DSM/IRFM) (yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM) (joan.decker@cea.fr)
0023 %
0024 if nargin < 2,
0025     error('Wrong number of input arguments in randreject_yp.m');
0026 elseif nargin < 3,
0027     N = 100000;
0028     display_mode = 0;
0029 elseif nargin < 4,
0030     display_mode = 0;
0031 end
0032 %
0033 y = NaN(N,1);
0034 %
0035 if nargin == 4,
0036     if display_mode,
0037         if round(N/display_mode) < 100,
0038             disp(['WARNING: reduce the number of bins or increase the number of trials for a correct histogramm -> display_mode parameter set to zero.']);
0039             display_mode = 0;            
0040         end
0041     end
0042 end
0043 %
0044 arg_range = args_in{length(args_in)};
0045 %
0046 if isnan(arg_range) | isempty(arg_range), 
0047     error('Error: the range of the input argument must be given');
0048 end
0049 %
0050 arg_index = length(args_in)+1;
0051 %
0052 if arg_range(1) < arg_range(2),
0053     xmax = fminbnd(@(x) (-1)*fun_handle(args_in(1:length(args_in)-1),[arg_range(1),arg_range(2)],x),arg_range(1),arg_range(2));
0054 else
0055     xmax = fminbnd(@(x) (-1)*fun_handle(args_in(1:length(args_in)-1),[arg_range(1),arg_range(2)],x),arg_range(2),arg_range(1));
0056 end  
0057 %
0058 if ~isempty(xmax)
0059     c = feval(fun_handle,args_in(1:length(args_in)-1),[arg_range(1),arg_range(2)],xmax); 
0060 else
0061     error('Error in determining the maximum of the function in the interval in randreject_yp.m.')
0062 end
0063 %
0064 while isnan(sum(y)),
0065     %
0066     y0 = (arg_range(2) - arg_range(1))*NaN(sum(isnan(y)),1) + arg_range(1);
0067     y1 = (arg_range(2) - arg_range(1))*rand(sum(isnan(y)),1) + arg_range(1);
0068     y2 = rand(sum(isnan(y)),1);
0069     %
0070     mask0 = y2 < feval(fun_handle,args_in(1:length(args_in)-1),[arg_range(1),arg_range(2)],y1)/c; 
0071     %
0072     y0(mask0) = y1(mask0);
0073     %
0074     inan = find(isnan(y) == 1);
0075     y(inan) = y0;
0076 end
0077 %
0078 if display_mode,
0079     y_th = (arg_range(2) - arg_range(1))*linspace(0,1,N) + arg_range(1);
0080     py_th = feval(fun_handle,args_in(1:length(args_in)-1),[arg_range(1),arg_range(2)],y_th);
0081     norm_th = trapz(y_th,py_th);
0082     %
0083     if round(N/display_mode) >= 100,%for good statistics, display_mode is the number of bin for the histogram
0084         [h1,y1] = hist(y,display_mode);
0085         py = h1*display_mode/N/(max(y1) - min(y1));
0086         norm = trapz(y1,py);
0087         close all
0088         figure(1),plot(y1,py,'b',y_th,py_th,'r')
0089         axis([min(arg_range(1),arg_range(2)),max(arg_range(1),arg_range(2)),min(py),max(py)])
0090         xlabel('y');ylabel('f(y)');
0091         if arg_range(2) > arg_range(1),
0092             title(['Distribution function between ',num2str(arg_range(1)),' and ',num2str(arg_range(2))])
0093             legend(['Numerical, norm. = ',num2str(norm)],['Theoretical, norm. = ',num2str(norm_th)],'Location','SouthWest')
0094         else
0095             title(['Distribution function between ',num2str(arg_range(1)),' and ',num2str(arg_range(2))])
0096             legend(['Numerical, norm. = ',num2str(norm)],['Theoretical, norm. = ',num2str(norm_th)],'Location','NorthWest')
0097         end
0098     else
0099         warning('Not enough points for histogram display');
0100     end
0101 else
0102     y_th = NaN;
0103     py_th = NaN;
0104     norm_th = NaN;
0105 end