Для решения поставленной я использовал среду программирования Delphi 7.
Сначала были сделаны расчеты по формулам которые описывались выше для того, чтобы отобразить график плотности логнормального распределения при помощи аналитических расчетов.
Для того, чтобы получить график плотности распределения на основе стохастических преобразований, был выбран метод Неймана.
Также нужно было посчитать математическое ожидание и дисперсию и ввести полученные результаты на экран. Формулы которых представлены ниже
| Математическое ожидание | 
|
| Дисперсия | 
|
var
Form1: TForm1;
kk:Int64;
flag:boolean;
implementation
Плотность распределения
function TForm1.PL(x:double):double; //--density of distribution
begin
if x<>0 then
result:= exp(-(ln(x)-mu)*(ln(x)-mu)/(2*sigma*sigma))/(x*sigma*Sqrt(2*Pi))
else
result:= 0;
end;
function TForm1.LogNorm(): double; //--for a method of Neumann
var
y: real;
x: double;
begin
repeat
x:= a+random*(b-a);
f:= PL(x);
y:= fmax*random;
until y<f;
result:= x;
end;
procedure TForm1.Clear; //------------clear array---------
const M=50;
var j: integer;
begin
for j:=0 to (M-1) do
begin
gist[j]:= 0;
end;
end;
procedure TForm1.Panel1Click(Sender: TObject);
var
x, r, sr, h1, h2, Ob,g1,g2, chi2_N, chi2_12, chi2_if, sum, Z: double;
p, y, Mat, Mat2, Disp: real;
M, j: integer;
N, i, u: longint;
begin
flag:=false;
Gauge1.Progress:=0;
//-------**All fields must be filled!**---------
if (E1.Text='') or (E2.Text='') or (E3.Text='') or (E4.Text='') or
(E5.Text='') then
begin
with Application do
begin
NormalizeTopMosts;
MessageBox('All of fields must be filled!', 'Error', MB_OK);
RestoreTopMosts;
end;
exit;
end;
//----------**initialization**--------------
T:= GetTime;
Clear;
Chart1.Series[0].Clear;
Chart1.Series[1].Clear;
Chart1.Series[2].Clear;
sigma:= StrToFloat(E1.Text);
mu:= StrToFloat(E2.Text);
a:= StrToFloat(E3.Text);
b:= StrToFloat(E4.Text);
kk:=StrToint64(E5.Text);
if kk>2000000000 then
begin
Showmessage ('Очень большое число, введите меньшее');
exit;
end;
N:= StrToInt(E5.Text);
g1:=100/N;
g2:=0;
Randomize;
M:= 50;
//---------------**theoretical method**------------------
for i:=1 to 100 do
begin
if (i mod 10) =0 then application.ProcessMessages;
x:= a+i*(b-a)/100;
//p:= PL(x);
if x<>0 then
p:= exp(-(ln(x)-mu)*(ln(x)-mu)/(2*sigma*sigma))/(x*sigma*Sqrt(2*Pi))
else
p:= 0;
Chart1.Series[0].AddXY(x, p);
end; //---theoretical
//***********************************************************
Метод Неймана
//---------------**method of Neumann**------------------------
fmax:=Chart1.Series[0].MaxYValue;
{for i:=1 to N do
begin
if (i mod 10) =0 then application.ProcessMessages;
x:= a+i*(b-a)/N;
f:= PL(x);
if (f>fmax)then
fmax:= f;
end;} //max
//------------------------------
Clear;
chi2_N:=0;
Mat:=0;
Mat2:=0;
Disp:=0;
i:=0;
Clear;
chi2_if:= 0;
while true do
begin
if (i mod 10) =0 then application.ProcessMessages;
inc(i);
x:= LogNorm();
Mat:= Mat+x; //expectation
Mat2:= Mat2 +sqr(x);
if (x>b) or (x<a) then
continue;
u:= trunc((x-a)/((b-a)/M));
gist[u]:= gist[u]+1;
h1:= random;
h2:= random;
Ob:= sqrt(-2*ln(h1))*cos(2*Pi*h2);
Ob:= mu+Ob*sigma;
x:= exp(Ob);
if (x>b) or (x<a) then
continue;
u:= trunc((x-a)/((b-a)/M));
gist1[u]:= gist1[u]+1;
g2:=g2+g1;
Gauge1.Progress:=trunc(g2)+1;
if i>N then break;
if flag=true then
begin
N:=i;
break;
end;
end;
Mat:= Mat/N;
Mat2:= Mat2/N;
Disp:= Mat2 - sqr(Mat);
for j:=0 to (M-1) do //------histogram
begin
sum:= (Power(N*PL(a+(b-a)/M*(j+0.5))*(b-a)/M-gist[j], 2))/
(N*PL(a+(b-a)/M*(j+0.5))*(b-a)/M);
chi2_N:= chi2_N+sum;
Chart1.Series[1].AddXY((a+(j+0.5)*(b-a)/M), gist[j]/N*M/(b-a));
end;
E6.Text:= FloatToStrF(chi2_N, fffixed, 4, 4);//--chi-square for a Neumann
//****************************************************************
Метод обратной функции
//--------------**method of inverse function**-----------------
Clear;
chi2_if:= 0;
{for i:=1 to N do
begin
h1:= random;
h2:= random;
Ob:= sqrt(-2*ln(h1))*cos(2*Pi*h2);
Ob:= mu+Ob*sigma;
x:= exp(Ob);
if (x>b) or (x<a) then
continue;
u:= trunc((x-a)/((b-a)/M));
gist[u]:= gist[u]+1;
end;}
for j:=0 to (M-1) do //------histogram
begin
sum:= (Power(N*PL(a+(b-a)/M*(j+0.5))*(b-a)/M-gist1[j], 2))/
(N*PL(a+(b-a)/M*(j+0.5))*(b-a)/M);
chi2_if:= chi2_if+sum;
Chart1.Series[2].AddXY((a+(j+0.5)*(b-a)/M), gist1[j]/N*M/(b-a));
gist1[j]:=0;
end;
E8.Text:= FloatToStrF(chi2_if, fffixed, 4, 4); //chi-sq for a inverse function
E10.Text:= FloatToStr(exp(mu+sqr(sigma)/2)); //--expectation (teor)
E11.Text:= FloatToStr(Mat); //--expectation (experim)
E12.Text:= FloatToStr((exp(sqr(sigma))-1)*exp(2*mu+sqr(sigma)));
E13.Text:= FloatToStr(Disp);
D:= GetTime;
Z:= MilliSecondSpan(D, T);
e5.Text:=IntTostr(N);
Edit1.Text:= FloatToStrF(Z, fffixed, 6, 6);
//*****************************************************************
end;
procedure TForm1.Panel7Click(Sender: TObject);
begin
Close;
end;
procedure TForm1.E1KeyPress(Sender: TObject; var Key: Char);
begin
if (key='-')
then begin
if Pos ('-', (Sender as TEdit).Text)=0 then Begin (Sender as TEdit).SelStart:=0; key:='-'; end
else key:=#0;
end;
if Sender is TEdit then
begin
if Not((Key in ['0'..'9'])or (Key=Chr(vk_Back))
or (Key=DecimalSeparator) or (Key='-')) then
Key:=#0
else
begin
if Key = DecimalSeparator then
if Pos(DecimalSeparator,(Sender as TEdit).Text)>0 then
Key:=#0;
end;
end;
end;
procedure TForm1.Aboutme1Click(Sender: TObject);
begin
AboutBox.Show;
end;
procedure TForm1.Timer1Timer(Sender: TObject);
begin
Panel19.Caption:= TimeToStr(Time);
end;
procedure TForm1.E1KeyDown(Sender: TObject; var Key: Word;
Shift: TShiftState);
begin
if (ssShift in Shift)then
key:=0;
end;
procedure TForm1.Panel20Click(Sender: TObject);
begin
flag:=true;
end;
end.
Инструкция пользователю
При запуске программы перед пользователем открывается форма, на которой есть поля ввода параметров, поля вывода посчитанных значений, поле для вывода графика и кнопки, при нажатии на которые происходит то или иное событие.
Справа в разделе "Теоретически пользователь может ввести значение sigma и mu, те значения которые он считает нужными; a и b это интервал в пределах которого меняется график. И значение N – (количество єксперементов) – в зависемости от того сколько раз мы будем проводить єксперемент. В зависимости от выбора данных параметров пользователь может получить различные формы графика плотности вероятности.
В разделе "Критерий согласия" выводятся значения оценки Хи-квадрат для двух указанных методов. Ниже вывод математического ожидания и дисперсии, посчитанных теоретически и экспериментально.
Справа внизу формы выводится системное время и время выполнения расчётов в миллисекундах.
При нажатии на кнопку «Вывести графики и вычислить» слева выводятся график плотности логнормального распределения (построенный теоретически), гистограммы распределения случайной величины по логнормальному закону, смоделированные при помощи метода Неймана и метода обратной функции.
При нажатии на кнопку «Стоп» программа прекращаются свою работу и начинает считывать значения которые обработались до определенного момента и записует значения в поля.
При нажатии на кнопку «2D/3D» пользователь может наблюдать изменение графика из 2D в 3D и наоборот.
При нажатии на кнопку «Выход» программа будет завершена.
В закладке «About» пользователь может узнать о создателях данного программного продукта и краткое описание программного продукта.