Текст подпрограммы и версий
gsr2r_p.zip 		Тексты тестовых примеров
tgsr2r_p.zip

Подпрограмма:  GSR2R (модуль GSR2R_p)

Назначение

Генерация массива псевдослучайных чисел, имеющих произвольное заданное распределение.

Математическое описание

По таблично заданной желаемой функции распределения F (x) строится интерполяционный многочлен  p (x) для функции F - 1 (x) - обратной по отношению к F (x).

Псевдослучайные величины с функцией распределения F (x) получаются по формуле R = p (u), где u - псевдослучайная величина, pавномеpно распределенная в интервале (0, 1).

Akima H., "A new method of interpolation and smooth curve fitting based on local procedures", JACM, 17(4), 1970, 589 - 602.

Guerra V.M., Tapia R.A., and Thomson J.R., "A random number generator for continuous random variables", ICSA Technical Report, Rice Univercity, Houston, Texas.

Использование 

procedure GSR2R(var ISEED :Integer; N :Integer; var R :Array of Real;
                var TBL :Array of Real; NX :Integer; var IP :Integer;
                var IERR :Integer);


   Параметры
ISEED - целая переменная, значение которой перед обращением к подпрограмме должно быть в пределах [1, 2147483646]; по окончании работы подпрограммы ей присваивается новое значение, котоpое может быть использовано при повторном вхождении в подпрограмму;
N - заданное количество псевдослучайных чисел, которые необходимо сгенерировать (тип: целый);
R - вещественный вектоp длины N содержащий полученные псевдослучайные величины;
T - двумерный вещественный массив размерностью NX * 5, содержащий табличные значения функции распределения F: элементы T (1, 1), T (2, 1), ..., T (NX, 1) содержат абсциссы, элементы T (1, 2), ..., T (NX,2) - ординаты (т.е. F (T (I, 1)) = T (I, 2)). Должны выполняться следующие условия: T (I, 1) < T (I + 1, 1), T (I, 2) < T (I + 1, 2) для I = 1, 2, ...,NX - 1; T (1, 2) = 0.; T (NX, 2) = 1.; остальные элементы массива используются для хранения промежуточных результатов;
NX - заданное количество узлов сетки, на которой задана функция распределения, NX ≥ 4 (тип: целый);
IP - задает режим работы подпрограмы (тип: целый); при этом если:
IP= 0 - то считается, что происходит первое обращение к подпрограмме и вычисляются коэффициенты интерполяционных полиномов;
IP= 1 - то коэффициенты интерполяционных полиномов считаются уже вычисленными;
  по окончании работы подпрограммы параметру IP присваивается значение единица;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR=65 - когда заданное количество узлов сетки, на которой задана функция распределения, меньше четырех;
IERR=66 - когда либо абсциссы, либо ординаты табличных значений функции распределения не являются строго монотонно возрастающими;
IERR=67 - когда T (1, 2) > 0, либо T (NX, 2) < 1. 
   Версии:  нет 
   Вызываемые подпрограммы
GSU1R - генерация массива псевдослучайных чисел, pавномеpно распределенных в интервале (0, 1);
I IS9R - интерполяция кубическим квази - эрмитовым сплайном табличной функции от одной переменной на неравномерной сетке;
UTGS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы GSR2R.
UTI I10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы I IS9R. 
   Замечания по использованию
  Если после обращения к GSR2R с одной функцией распределения необходимо обратиться к GSR2R с другой функцией распределения, то необходимо установить параметр IP = 0.

Пример использования 

Unit TGSR2R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, GSR2R_p;

function TGSR2R: String; 

implementation

function TGSR2R: String;
var
NX,N,IP,ISEED,_i,IERR :Integer;
R :АRray [0..2] of Real;
const
T :Array [0..29] of Real = ( 0.0,0.6,2.0,4.2,7.2,11.0,0.0,0.2,0.4,0.6,0.8,
1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0 );
begin
ReSULt := '';
NX := 6;
N := 3;
IP := 0;
ISEED := 12457;
GSR2R(ISEED,N,R,T,NX,IP,IERR);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[R[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%5d %5d',[IP,IERR]) + #$0D#$0A;
UtRes('TGSR2R',Result);  { вывод результатов в файл TGSR2R.res }
exit;
end;

end.

Результаты:

       IP    =  1
       R     =  ( 0.192542,  3.772356,  3.887485 )
       IERR  =  0