randlin_yp

 Calculate the random variable with a linear decrease of the p.d.f. (normalized to 1 over the the domain).

 Input:

   - N : number of trials [1,1]
     (default = 100000)
   - h : ratio between the pdf values at ybounds [1,1]
     (default = 1)
   - ybounds : range of the random variable [3,1] or [1,3] for example. If [1,3], f(1)/f(3) = 1/h, and if [3,1], f(1)/f(3) = h. If NaN or empty ybounds = [0,1]
     (default [0,1])
   - 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] = randlin_yp(N,h,ybounds,display_mode)
0002 %
0003 % Calculate the random variable with a linear decrease of the p.d.f. (normalized to 1 over the the domain).
0004 %
0005 % Input:
0006 %
0007 %   - N : number of trials [1,1]
0008 %     (default = 100000)
0009 %   - h : ratio between the pdf values at ybounds [1,1]
0010 %     (default = 1)
0011 %   - ybounds : range of the random variable [3,1] or [1,3] for example. If [1,3], f(1)/f(3) = 1/h, and if [3,1], f(1)/f(3) = h. If NaN or empty ybounds = [0,1]
0012 %     (default [0,1])
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 < 1, 
0025     N = 100000;
0026     h = 1;
0027     ybounds = [0,1];
0028     display_mode = 0;
0029 elseif nargin < 2, 
0030     h = 1;
0031     ybounds = [0,1];
0032     display_mode = 0;    
0033 elseif nargin < 3,
0034     ybounds = [0,1];
0035     display_mode = 0;
0036 elseif nargin < 4,
0037     display_mode = 0;
0038 end
0039 %
0040 if isnan(ybounds) | isempty(ybounds), 
0041     ybounds = [0,1];
0042 end
0043 %
0044 if h == 1,%uniform distribution with p(0) = 1 and p(1) = 1
0045     x = rand(N,1);
0046 elseif h == 0,%linear decreasing distribution with p(0) = 1 and p(1) = 0
0047     x1 = rand(N,1);
0048     x2 = rand(N,1);
0049     %
0050     x = abs(x2-x1);
0051 else,%linear decreasing distribution with p(0) = 1 and p(1) = h (rejection technique)
0052     %
0053     x = NaN(N,1);
0054     %
0055     while isnan(sum(x)),
0056         x1 = rand(sum(isnan(x)),1);
0057         x2 = rand(sum(isnan(x)),1);
0058         x0 = NaN(sum(isnan(x)),1);
0059         %
0060         mask0 = x2 < (1 - (1-h)*x1);
0061         %
0062         x0(mask0) = x1(mask0);
0063         %
0064         inan = find(isnan(x) == 1);
0065         x(inan) = x0;
0066     end
0067 end
0068 %
0069 if display_mode < 10;
0070     display_mode = 10;
0071     warning('display_mode must exceed 10 for a correct histogramm.');
0072 end
0073 %
0074 y = (ybounds(2) - ybounds(1))*x + ybounds(1);
0075 %
0076 y_th = sort((ybounds(2) - ybounds(1))*linspace(0,1,N) + ybounds(1));%increasing order
0077 py_th = (1 - (1 - h)*(y_th - ybounds(1))/(ybounds(2) - ybounds(1)))/abs((ybounds(2) - ybounds(1)))*2/(h+1);%p(y) = p(x)*abs(dx/dy)
0078 norm_th = trapz(y_th,py_th);
0079 %
0080 if display_mode,
0081     if N/display_mode > display_mode,%for good statistics, display_mode is the number of bin for the histogram
0082         [h1,y1] = hist(y,display_mode);
0083         py = h1*display_mode/N/(max(y1) - min(y1));
0084         norm = trapz(y1,py);
0085         close all
0086         figure(1),plot(y1,py,'b',y_th,py_th,'r')
0087         axis([min(ybounds(1),ybounds(2)),max(ybounds(1),ybounds(2)),min(py_th),max(py_th)])
0088         xlabel('y');ylabel('f(y)');
0089         if ybounds(2) > ybounds(1),
0090             title(['Linear distribution function between ',num2str(ybounds(1)),' and ',num2str(ybounds(2)),' -> f(',num2str(ybounds(1)),')/f(',num2str(ybounds(2)),') = ',num2str(1/h)])
0091             legend(['Numerical, norm. = ',num2str(norm)],['Theoretical, norm. = ',num2str(norm_th)],'Location','SouthWest')
0092         else
0093             title(['Linear distribution function between ',num2str(ybounds(1)),' and ',num2str(ybounds(2)),' -> f(',num2str(ybounds(1)),')/f(',num2str(ybounds(2)),') = ',num2str(h)])
0094             legend(['Numerical, norm. = ',num2str(norm)],['Theoretical, norm. = ',num2str(norm_th)],'Location','NorthWest')
0095         end
0096     else
0097         warning('Not enough points for histogram display -> set display_mode to ');
0098     end
0099 end