lw2

clear;
clc;
%{
1. В соответствии с номером варианта, используя генератор псевдослучайных чисел
(встроенную в математический пакет функцию),
генерировать выборку из заданного непрерывного распределения объемом N=10000 с заданными для вашего распределения
параметрами
нормальный закон с параметрами µ=1 σ=2
%}

imu = 1;
isigma = 2;
N = 1000;

type = "Normal";
X = random(type, imu, isigma, N, 1);

%{
2. Получение точечных оценок по заданной выборочной совокупности.
2.1. По заданной совокупности выборочных значений найти оценку параметра заданного распределения,
из которого генерировалась выборка, используя метод моментов. Для этого:
-составить уравнение, приравняв какой-либо теоретический момент соответствующему эмпирическому моменту того же порядка
- решить получившиеся уравнение относительного оцениваемого вами параметра
-проанализировать поведение оценки параметра при изменении объема выборки
-сделать выводы
%}

syms cmu csigma;

%функции для решения в уравнениях (приравниваем теоретический момент эмпирическому соответствующего порядка)
fun_mu = cmu == mean(X);
fun_sigma = csigma^2 == moment(X, 2);

%оценочное мю и сигма
ans_mu = double(solve(fun_mu));
ans_sigma = double(max(solve(fun_sigma)));

%{
2.2. По заданной совокупности выборочных значений найти оценку параметра заданного распределения,
из которого генерировалась выборка, используя метод максимального правдоподобия. Для этого необходимо:
- построить график функции правдоподобия от оцениваемого параметра
- найти положение максимума этой функции. Это значение и будет точечной оценкой
- проанализировать поведение оценки при изменении объема выборки
-проанализировать поведение функции правдоподобия при изменении объема выборки
-сделать выводы
%}

syms sigma;

%этот код работает только для маленьких чисел и является более честной
%версией генерации следующей функции

% por = sym('a', [1 N]);
%
% for i = 1:N
% por(i) = (1 / (sigma * sqrt(2 * pi))) * exp ( -(X(i) - imu)^2 / ( 2 * sigma^2 ) );
% end
% clear i;
%
% pfun = prod(por);

%создаем функцию
pfun = (1 / (sigma * sqrt(2 * pi)))^N * exp ( sum( -(X - imu) .^2 ./ ( 2 * sigma^2 ) ) );

%ищем логарифм от функции
Lg= log(pfun);
%получаем значения, в которых производная равна нулю ( максимум )
dLg = diff(Lg, sigma);
eq0_dLg = dLg == 0;

%построение функции, работает только для N ~ 1000
figure;
hold on;
fplot(Lg, sigma);
hold off;

%производим оценку параметра
ans_solve_sigma = abs(double(vpasolve(dLg, sigma)));

%{
2.3. Оценить дисперсию, смещение и рассеяние оценки параметра по 100 выборочным реализациям
%}

n = 100;

mean_100 = linspace(1, n);
var_100 = linspace(1, n);

for i = 1:n

    rng shuffle;
    X_100 = random(type, imu, isigma, N, 1);

    mean_100(i) = 1/mean(X_100);
    var_100(i) = moment(X_100, 2);

end
clear i;

ans_mean_m = mean(mean_100) - mean(X);
ans_mean_v = mean(var_100) - mean(X);
ans_var_m = var(mean_100);
ans_var_v = var(var_100);
ans_rass_m = mean((mean_100 - mean(X)).^2);
ans_rass_v = mean((var_100 - mean(X)).^2);

clear n;

%генерация двух дополнительных выборок с использованием параметров,
%полученных в результате точечных оцениваний:

X_sigmacustom1 = random(type, imu, ans_sigma, 1, N);
X_sigmacustom2 = random(type, imu, ans_solve_sigma, 1, N);

%{
3. Расчет гистограммы относительных частот.
3.1. Построить гистограмму по уже описанному в первой лабораторной работе алгоритму
%}

X_max = max(X);
X_min = min(X);

r = floor(log2(N) + 1);
h = (X_max - X_min) / r;

z = zeros(0,r);

for i = 0:r
    z(i+1) = X_min + i*h;
end

clear i;

zX = zeros(1, r);

for i = 1:r
    zX(i) = z(i+1) - h/2;
end

clear i z;

U = hist(X, zX);

H_X = U / (h*N);

%{
3.2. Нанести на один и тот же график плотность вероятности того теоретического распределения,
из которого генерируется выборка
3.3 Нанести на тот же график плотность вероятности того теоретического распределения,
из которого генерируется выборка но вместо истинного значения параметра, подставить оценку,
полученную ранее с использованием одного из методов точечного оценивания.
%}

max_X = max(X);

x = 0:0.1:max_X;

y = pdf(type, x, imu, isigma);

y1 = pdf(type, x, imu, ans_sigma);

y2 = pdf(type, x, imu, ans_solve_sigma);

figure
hold on
bar(zX, H_X, 1);
plot(x,y,':','Color',[0 0.4470 0.7410],'LineWidth',1.5);
plot(x, y1,'--.','Color',[1 0 1],'LineWidth',1.5);
plot(x, y2,'-.','Color',[0.6350 0.0780 0.1840],'LineWidth',1.5);
hold off
title('Пункты 3.2, 3.3');
legend({'Гистограмма', 'Плотность', 'Плотность(сигма 1)', 'Плотность(сигма 2)',}, 'Location','northeast');
grid on
%{
3.4 Сделать выводы
%}

%{
4. Расчет эмпирической функции распределения.
Выполнить расчет эмпирической функции распределения для негруппированных данных.
Также, необходимо рассчитать и построить вместе с этой функцией график теоретической функции того распределения,
из которого выборка генерировалась с заданными в условии параметрами.
И построить график теоретической функции распределения, подставив вместо истинного значения параметра его оценку.
Сравнить все представленные функции и сделать выводы
%}

y_cdf1 = cdf(type, x, imu, isigma);
y_cdf2 = cdf(type, x, imu, ans_sigma);
y_cdf3 = cdf(type, x, imu, ans_solve_sigma);

figure;
hold on
plot(x,y_cdf1,':','Color',[0 0.4470 0.7410],'LineWidth',1.5);
plot(x, y_cdf2,'--.','Color',[1 0 1],'LineWidth',1.5);
plot(x, y_cdf3,'-.','Color',[0.6350 0.0780 0.1840],'LineWidth',1.5);
hold off
title('Эмпирическая функция распределения');
legend({'Стандартная', 'Измененная сигма 1', 'Измененная сигма 2'},'Location','southwest')