Генерация отсчетов случайной величины с заданной плотностью распределения методом кусочно-линейной аппроксимации

 % msv_kla.m
% Раздел 3. 
% Моделирование С.В. методом кусочно-линейной апроксимации плотности
% распределения (при невозможности отыскания функции, обратной закону распределения)
%----------------------------------------------------------------------
clc, clear
 
kk=10000;   %input(' :> kk=');  % размер генеральной выборки
n=100;                          % число интервалов для отображения на гистограмме
 
% выбор варианта исходных плотностей ----------------------------------
prForm=1; % 1 - gaus; 2 - reley; 3 - <file>
switch prForm
    case 1, %формирование плотности распределения в программе
        Mx=10; Dx=1; % задаем параметры нормального закона распределения
        xmin=Mx-3*sqrt(Dx); xmax=Mx+3*sqrt(Dx); dx=(xmax-xmin)/n;
        xmax=xmin+n*dx;  %коррекция максим.границы
        x=xmin:dx:xmax;
        fx=exp(-((x-Mx).^2/(2*Dx)))/sqrt(2*pi*Dx);
    case 2, %формирование плотности распределения в программе
        % задаем параметры закона Рэлея (для разнообразия)
        sigm=3; %параметр масштаба или мода плотности распределения
        xmin=0;     xmax=11; %11->sigm=3; 15->sigm=4; 8->sigm=2
        dx=(xmax-xmin)/n;
        xmax=xmin+n*dx;  %коррекция максим.границы
        x=xmin:dx:xmax;
        for i=1:length(x), 
           fx(i)=x(i)/sigm^2*exp(-x(i)^2/(2*sigm^2));
        end;
    case 3, %загрузка плотности распределения из файла (подготавливается заранее)
        load file_fx;   % в файле x,fx,dx, !!! kk, n
    otherwise disp('данные не подготовлены'); return
end;
 
figure(1);plot(x,fx);grid;  % график "точной" плотности распределения
 
%-------------------------------------------------------------------
% Моделирование С.В. методом кусочно-линейной апроксимации плотности
%-------------------------------------------------------------------
% определение вероятностей попадания в полученные интервалы
for i=1:n,
    fxsr(i)=(fx(i+1)+fx(i))/2; p(i)=fxsr(i)*dx;
end;
% нормирование (для провед единичного испытания)
s=zeros(1,n); sp=0; sumP=sum(p);
for i=1:n, 
    pn(i)=p(i)/sumP; 
    sp=sp+pn(i); spx(i)=sp;
end;
% собственно моделирование СВ, равномерно распределенной на интервале [x(i)...x(i+1)[
for ik=1:kk,
    r1=rand;    % одиночное испытание - попадание в один из интервалов
    iI=1;    while r1>spx(iI), iI=iI+1; end;
    r2=rand; xm(ik)=x(iI)+r2*(x(iI+1)-x(iI));   % генерация равном.распр.СВ на интевале
end;
% визуальный контроль сформированной выборки по гистограмме ---------
figure(2);hist(xm,n);grid;  % гистограмма полученного распределения