Планирование и обработка результатов пассивного эксперимента

Теоретическая часть

Основная задача пассивного эксперимента — по результатам наблюдений сделать некоторые выводы о параметрах математической модели эксперимента [1]. При этом вид ее предполагается известным, а параметры — неизвестными. Далее будет рассматриваться класс линейных регрессионных моделей эксперимента.

В общем случае объект исследования представляется в виде схемы, представленной на рис. 9.1, где $$X = \{x_{1}, x_{2}, … , x_{k}\}$$ — входные величины или факторы; $$y_{i}$$ — $$i$$ -я выходная величина ( $$i = 1, 2, … , n$$ ); $$\epsilon$$ — случайные неконтролируемые возмущения [12].

Под моделью объекта по $$i$$ -му каналу понимают функцию

$$y_i = f(x_1, x_2, … , x_k).$$

Так как имеются случайные неконтролируемые возмущения, изменение функции (9.1) носит случайный характер, а потому для получения математического описания (9.1) применяются методы регрессионного анализа на основе статистических данных, накопленных в результате проведения эксперимента.

(рис 9.1) Схема объекта исследования

Применение метода пассивного эксперимента может быть успешным, если при его проведении соблюдаются необходимые условия, к которым относятся такие, как правильное определение времени регистрации данных, обеспечение независимости соседних измерений и входных переменных друг от друга, достаточный с точки зрения математической статистики объем экспериментальных данных [12].

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

$$y=f(x_1,x_2,...,x_k)=b_0+\sum\limits_{j=1}^{k}b_jx_j+\sum\limits_{i,j=1}^{k}b_{ij}x_ix_j+\sum\limits_{j=1}^{k}b_{jj}x_j^2+...,$$

где $$b_{0}$$, $$b_{j}$$, $$b_{ij}$$, $$b_{jj}$$ — постоянные коэффициенты уравнения, оценки которых необходимо определить в результате постановки и проведения пассивного эксперимента; $$k$$ — число наиболее существенных входных величин, полученных в результате отсеивающего эксперимента.

Пространство, в котором строится поверхность отклика – реакция выходной величины, называется факторным пространством .

Для применения методов регрессионного анализа требуется выполнение ряда предпосылок [12], а именно:

результаты наблюдений $$y_{1}, y_{2}, … , y_{N}$$ выходной величины в точках факторного пространства представляют собой независимые случайные величины, распределенные по нормальному закону, а процесс изменения выходной величины должен быть стационарным во времени;

дисперсии $$D_{yi} (i = 1, 2, … , N)$$ этих случайных величин должны быть равны друг другу (выборочные оценки дисперсий однородны);

все значения входных величин должны измеряться с пренебрежимо малой ошибкой по сравнению с ошибкой измерения выходной величины;

входные величины не должны коррелировать между собой;

все соседние измерения по каждой $$j$$ -й входной величине должны быть независимы.

Число коэффициентов уравнения (9.2) определяет объем эксперимента. Поэтому выбирают такой полином, который содержит как можно меньше коэффициентов, но удовлетворяет требованию простоты и адекватности, под которой понимается способность модели предсказывать результаты эксперимента в некоторой области с требуемой точностью.

Часто на предварительной стадии исследования объекта выбирают полином первой степени (включая первую сумму в уравнении (9.2), т. е. в разложении Тейлора), предполагая, что параметры объекта лежат в области, в которой расположен экстремум исследуемой функции (выходная величина системы – отклик), и поэтому объект описывается линейной моделью. Если же эта линейная модель оказывается неадекватной, то в нее включают члены парного взаимодействия $$x_{i}x_{j}$$ (включая вторую сумму в уравнении (9.2), т. е. в разложении Тейлора), а при необходимости увеличивают степень полинома до тех пор, пока модель не окажется адекватной [12].

В результате регрессионного анализа результатов пассивного эксперимента находят оценки коэффициентов уравнения регрессии $$\beta_{0},\mbox{ } \beta_{j},\mbox{ }\beta_{ij},\mbox{ }\beta_{jj},\mbox{ }…$$.

1. План эксперимента

Рассмотрим эксперимент, в котором проводится $$N$$ измерений зависимой переменной $$y$$ в некоторых точках факторного пространства [1]. Обозначим через $$y_{u}$$ наблюдаемое значение зависимой переменной $$у$$ в $$u$$ -м опыте в точке

$$х^{u} = (x_{1u}, x_{2u}, …, x_{ku}), u = 1, 2, … , N$$.

Здесь $$x_{iu} (i = 1, 2, … , k)$$ — значение переменной $$x_{i}$$ в $$u$$ -м опыте.

Определение 1. Набор точек $$x^{u}\mbox{ }(u = 1, 2, … , N)$$ называется планом эксперимента. Точки при этом не обязательно должны быть различными. Матрица вида

$$X=\left[\begin{array}{cccc} x_{11} x_{21} \ldots x_{k1}\\ x_{12} x_{22} \ldots x_{k2}\\ \ldots \ldots \ldots \ldots\\ x_{1N} x_{2N} \ldots x_{kN} \end{array} \right]$$

называется матрицей плана эксперимента.

Если обозначить различные точки плана через $$х_{1}, х_{2}, … , х_{n}$$, то совокупность таких точек называется спектром плана.

Определение 2. Нормированным планом $$\varepsilon(N)$$ называют совокупность величин $$х_{1}, х_{2}, … , х_{n}$$ ; $$р_{1}, р_{2}, … , р_{n}$$, где $$\sum\limits_{i=1}^{n}p_i,\mbox{ }p_i=m_i/N,\nbox{ }\sum\limits_{i=1}^{n}m_i=N$$.

В пассивном эксперименте задача построения плана не рассматривается. Матрица плана (9.3) предполагается известной (заданной) или является предопределенной условиями проведения эксперимента [1]. Задача исследователя в пассивном эксперименте состоит в выполнении наблюдений над выходной (зависимой) переменной в точках, определяемых матрицей плана, и последующем анализе их результатов.

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

$$X=\left[\begin{array}{ccccc} 1x_{11} x_{21} \ldots x_{k1}\\ 1x_{12} x_{22} \ldots x_{k2}\\ \ldots\ldots \ldots \ldots \ldots\\ 1x_{1N} x_{2N} \ldots x_{kN} \end{array} \right]$$

2. Одномерная регрессионная модель эксперимента

В общем случае вид функции отклика неизвестен. Будем предполагать, что функция отклика является одномерной и представима в виде

$$y=\sum\limits_{j=1}^{p}f_{j}(x_1,x_2,...,x_k)\beta_j,$$

где:

$$\beta = (\beta_1, \beta_2, … , \beta_р)$$ есть $$р$$ -мерный вектор неизвестных параметров;

$$\{f_{j}(x_{1}, x_{2}, … , x_{k})\}$$ — известные функции, которые называются еще базисными функциями.

Под регрессионной моделью эксперимента будем понимать линейную по параметрам $$\beta_1, \beta_2, … , \beta_р$$ функцию отклика (9.5).

В качестве примера рассмотрим функцию отклика

$$у = \beta_0 + \beta_1х + \beta_{11}х^2.$$

Построим матрицу планирования по заданным измерениям входной переменной $$х = (2, 4, –3, 1)$$. В соответствии с (9.4) будем иметь

$$X=\left[\begin{array}{ccccc} 12 4\\ 14 16\\ 1-3 9\\ 111 \end{array} \right] .$$

Еще один пример. Пусть функция отклика имеет вид

$$y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_{12}x_1x_2.$$

Значения независимых переменных равны следующим значениям:

$$х_{1} = (–1, 1, –1, 1),\mbox{ }х_{2} = (–1, –1, 1, 1)$$.

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

$$X=\left[\begin{array}{ccccc} x_0x_1 x_2 x_1x_2\\ 1-1 -1 1\\ 11 -1 -1\\ 1-1 1 -1\\ 11 1 1\\ \end{array} \right],$$

где $$х_{0} = 1$$ — фиктивная переменная.

3. Оценивание параметров одномерной функции отклика

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

$$y=\sum\limits_{j=1}^{p}f_{j}(x_1,x_2,...,x_k)\beta_j,$$

в произвольной точке $$(x_1, x_2, … , x_k)\in G\subset R^k$$ по результатам наблюдений (измерений) в точках $$x_{1}, x_{2}, … , x_{N}$$, где $$x_{u} = (x_{1u}, x_{2u}, … , x_{ku})$$, задаваемых матрицей плана (планирования) $$M = (x_{iu}) (i = 1, 2, … , k; u = 1, 2, … , N)$$, наилучшая оценка в смысле метода наименьших квадратов равна

$$\hat\beta=(X^TX)^{-1}X^TY,$$

если ранг матрицы планирования равен числу неизвестных параметров, т. е. $$rank(X) = p$$, где $$р$$ — число неизвестных параметров [1].

Оценка параметров по формуле (9.6) называется оценкой по наблюдениям полного ранга, т. е. когда $$rank(X) = p$$. В этом случае формула (9.6) вытекает из так называемого нормального уравнения

$$(X^TX)\hat\beta=X^TY,$$

где $$X^{T}X$$ — информационная матрица.

В методе наименьших квадратов минимизируется функционал (скалярная величина)

$$Q = (Y – Xb)^T(Y – Xb),$$

где $$b = (b_{1}, b_{2}, … , b_{p})$$ — вектор неизвестных параметров, подлежащих оценке.

Оценка параметров по формуле (9.6) доставляет минимум функционалу (9.8). При подстановке (9.6) в (9.8) получающаяся величина называется остаточной суммой квадратов — RSS (Residual Sum of Squares).

4. Оценивание параметров многомерной функции отклика

В случае, когда имеются несколько выходных переменных — функций отклика, для каждой функции отклика можно произвести оценку параметров уравнения регрессии по формуле (9.6) с указанием нижнего индекса, относящегося к данной функции. Например, если имеется $$n$$ функций отклика, то оценка вектора параметров по методу наименьших квадратов может быть найдена как

$$\hat\beta=(X^T_jX_j)^{-1}X^T_jY_j,\qquad j=1,2,...,n.$$

Общая оценка вектора параметров многомерной функции отклика будет определяться в виде

$$\hat\beta=(\hat\beta_1^T,\hat\beta_2^T,...,\hat\beta_n^T),$$

где $$Т$$ — символ транспонирования.

Практическая часть

1. Оценка параметров одномерной функции отклика

При создании моделирующего алгоритма оценивания параметров одномерной функции отклика принимаем во внимание, что наблюдения функции отклика представляют собой случайные величины, распределенные по нормальному закону [1, 12]. При этом ранг матрицы планирования будем считать равным числу оцениваемых параметров. Входные величины в принципе могут быть произвольными.

Моделирующий алгоритм будет следующим:

Задание базы выхода.

Задание входных воздействий по произвольному вероятностному закону.

Определение вектора выхода для r опытов.

Определение матрицы входа для r опытов.

Определение оценки вектора параметров.

Проверка заданной погрешности.

Расчет остаточной суммы квадратов.

Пример 1. По заданной базе выхода (например, 6.78) и трем входным воздействий, заданных по вероятностному закону, определите параметры линейной регрессионной зависимости. Произведите оценку погрешности и рассчитайте остаточную сумму квадратов.

Программный код решения примера:

clear,clc
% ПАРАМЕТРЫ ДИАЛОГОВОГО ОКНА inputdlg
options.Resize = 'on';
options.WindowStyle = 'normal';
options.Interpreter = 'tex';
%--------------------------------------
X = inputdlg({'\bf Ввод значения функции отклика...................'}, 'База функции отклика', 1, {'6.78'}, options);
Xb = str2num(char(X));
if prod(size(Xb)) ~= 1 | isempty(Xb) == 1 | ...
        ~isreal(Xb) | isinf(Xb) | isnan(Xb)
    errordlg('Отклик функции должен быть одним действительным числом',...
        'Ошибка ввода функции отклика');
    return
end
ms = inputdlg({'\bfМатематическое ожидание нормального закона:','\bfСтандартное отклонение нормального закона: .......'},...
    'Параметры',1,{'0','1'},options);
m = str2num(char(ms(1)));
s = str2num(char(ms(2)));
if isempty(m) | prod(size(m)) ~= 1 | ~isreal(m) | isinf(m) | isnan(m)
    errordlg('Математическое ожидание должно быть одним действительным числом',...
        'Ошибка ввода математического ожидания');
    return
end
if isempty(s) | prod(size(s)) ~= 1 | ~isreal(s) | ...
        isinf(s) == 1 | s <= 0 | isnan(s)
    errordlg('Стандартное отклонение должно быть действительным числом больше нуля',...
    'Ошибка ввода стандартного отклонения');
    return
end
% ВВОД РАЗМЕРНОСТИ ВХОДА И ЧИСЛА ОПЫТОВ
u = inputdlg({'\bfЧисло входов:...',...
    '\bfЧисло опытов в эксперименте:..................................'},...
 'Параметры эксперимента:',1,{'3','10'},options);
p = str2num(char(u));
if isempty(p) | prod(size(p(1))) ~= 1 |  ~isreal(p(1)) | ...
    isinf(p(1)) | p(1) ~= round(p(1)) | p(1) <= 0 | isnan(p(1))
    errordlg('Число входов должно быть натуральным числом',...
        'Ошибка ввода числа входов');
    return
end
if isempty(p(2)) | prod(size(p(2))) ~= 1 | isnan(p(2)) | ...
       ~isreal(p(2)) | isinf(p(2)) | p(2) ~= round(p(2)) 
    errordlg('Число опытов должно быть натуральным числом',...
        'Ошибка ввода числа опытов');
    return
end
if p(2) <= p(1) + 1
errordlg('Число опытов должно быть больше числа входов на 2 единицы ',...
    'Ошибка ввода числа опытов');
    return
end
%-----------------------------------------------------------------
% ЗАДАНИЕ ПОГРЕШНОСТИ ПО ВЫХОДУ
if abs(Xb) > 0.001
zet = inputdlg({'\bfВеличина относительной погрешности в процентах:..................'},...
 'Относительная погрешность по выходу',1,{'0.5'}, options);
dz = str2num(char(zet));
if isempty(dz) | prod(size(dz)) ~= 1 | ~isreal(dz) | ...
        isinf(dz) == 1 | dz <= 0 | isnan(dz)
    errordlg('Погрешность должна быть действительным числом больше нуля...',...
            'Ошибка ввода относительной погрешности');
    return
end

else
   zet = inputdlg({'\bfВеличина абсолютной погрешности по выходу:...'},...
       sprintf('Абсолютная погрешность '),1,{'0.001'}, options);
dz = str2num(char(zet));
if isempty(dz) | prod(size(dz)) ~= 1 | isreal(dz) == 0 | ...
        isinf(dz) == 1 | dz <= 0 | isnan(dz)
    errordlg('Погрешность должна быть действительным числом больше нуля',...
            'Ошибка ввода абсолютной погрешности');
    return
end
end
%-----------------------------------------------------------------
Ub = [1,randn(1,p(1))];
%------------------------------------------------------------------
k = p(2); % ЧИСЛО ОПЫТОВ
d = 100;
n = 0;
% dz - ПОГРЕШНОСТЬ ПО ВЫХОДУ
while d > dz
n = n + 1;
% ФОРМИРОВАНИЕ ВЕКТОРА ВЫХОДА С ПОВТОРНЫМИ ОПЫТАМИ
for I = 1:k
    R(I) = normrnd(m,s);
end
X = Xb + R';
%------------------------------------
U = Ub(2) + randn(k,1);
% p(1) - ЧИСЛО ВХОДОВ
for J = 1:p(1) - 1
    U = [U,Ub(J) + rand(k,1)];
end
U1 = [ones(length(U(:,1)),1),U];
%------------------------------------
b = inv(U1'*U1)*U1'*X;
%------------------------------------
% РАСЧЕТНОЕ ЗНАЧЕНИЕ ВЫХОДА
Xpac = Ub*b;
%------------------------------------
% ОТНОСИТЕЛЬНАЯ ПОГРЕШНОСТЬ В ПРОЦЕНТАХ
if abs(Xb) > 0.001

d = abs(Xpac-Xb)/abs(Xb)*100;
else
    d = abs(Xpac-Xb);
end
end
%------------------------------------
% ОСТАТОЧНАЯ СУММА КВАДРАТОВ
Q = (X - U1*b)'*(X - U1*b);
%------------------------------------
disp(sprintf('\n ОБРАБОТКА РЕЗУЛЬТАТОВ ПАССИВНОГО ЭКСПЕРИМЕНТА'));
disp(sprintf(' Базовый выход Xb = %g',Xb))  
fprintf(' Параметры нормального закона для формирования повторных опытов выхода:\n');
 
fprintf('\t математическое ожидание  m = %g\n',m);
fprintf('\t cтандартное отклонение  s = %g\n',s);
disp(sprintf('\t Число входов m = %g',length(Ub) - 1));
disp(sprintf('\t Значения входного воздействия:'));
 
for J = 2 : length(Ub)
    fprintf('\t %-g\n', Ub(J));
end
 fprintf('\n');
if abs(Xb) > 0.001
 fprintf('\t Заданная относительная погрешность в процентах: %g%s\n',dz,'%')
else
fprintf('\t Заданная абсолютная погрешность: %g\n',dz);
end
disp('--------------------------------------------------------------')
fprintf('\t Число повторных опытов в эксперименте: %d\n',p(2))
fprintf('\n\t Расчетная матрица входных воздействий с повторными опытами:\n');
disp(U1)
fprintf('\t Расчетный вектор выхода с повторными опытами:\n')
fprintf('\t %-g\n',X)
if abs(Xb) > 0.001
disp('--------------------------------------------------------------')
disp(sprintf('\t Расчетный выход Xpac = %g\n\t Погрешность расчета в процентах d = %g%s',...
    Xpac,d,'%'))
elseif abs(Xb) <= 0.001
disp('--------------------------------------------------------------')
disp(sprintf('\t Расчетный выход Xpac = %g\n\t Абсолютная погрешность расчета  d = %g', Xpac,d))
end
 fprintf('\t РАСЧЕТНЫЕ КОЭФФИЦИЕНТЫ УРАВНЕНИЯ РЕГРЕССИИ:');
for J = 1 : length(b)
    fprintf('\n\t b%d = %-g', J-1, b(J));
end
fprintf('\n');
 disp(sprintf('\t Остаточная сумма квадратов Q = %g',Q))
disp(sprintf('\t Число итераций по условию проверки погрешности n = %g',n))
fprintf('\n\t 1) Символьное уравнение регрессии:\n');
bn = length(b);
fprintf('\t\t y = '),fprintf('b0 + ');
for J = 1:bn - 2
    fprintf('b%d*x%d + ',J,J);
end
fprintf('b%d*x%d\n',bn - 1,bn - 1);
fprintf('\t 2) Уравнение регрессии с числовыми оценочными параметрами:\n')
bn = length(b);
fprintf(' y = '),fprintf('%g + ',b(1));
for J = 2 : (bn - 1)
    fprintf('(%g)*x%d + ',b(J),J - 1);
end
fprintf('(%g)*x%d\n',b(end), bn - 1);
disp('--------------------------------------------------------');

Ниже приводится одна из возможных реализаций выполнения программы (т. к. выход и входы смоделированы по случайным законам).

ОБРАБОТКА РЕЗУЛЬТАТОВ ПАССИВНОГО ЭКСПЕРИМЕНТА
 Базовый выход: Xb = 6.78
 Параметры нормального закона для формирования повторных опытов выхода:
	 математическое ожидание:  m = 0
	 cтандартное отклонение:  s = 1
	 Число входов: m = 3
	 Значения входного воздействия:
	 -1.01743
	 -1.4031
	 -0.260899

	 Заданная относительная погрешность в процентах: 0.5%
--------------------------------------------------------------
	 Число повторных опытов в эксперименте: 10
	 Расчетная матрица входных воздействий с повторными опытами:
    1.0000   -0.0259    1.8706   -0.8428
    1.0000   -2.0962    1.5201   -0.2298
    1.0000   -0.7237    1.4284   -0.1887
    1.0000   -1.4151    1.6469   -0.7250
    1.0000   -0.8569    1.3781   -0.1382
    1.0000   -0.3895    1.9245   -0.0297
    1.0000   -2.4837    1.0058   -0.9198
    1.0000    0.0256    1.7664   -0.6273
    1.0000   -1.1517    1.6920   -0.8416
    1.0000   -0.9158    1.1305   -0.5560
	 
        Расчетный вектор выхода с повторными опытами:
	 5.86991
	 5.30057
	 9.27119
	 9.05969
	 6.35003
	 9.75507
	 8.46856
	 6.92237
	 9.01667
	 5.847
--------------------------------------------------------------
	 Расчетный выход: Xpac = 6.80021
	 Погрешность расчета в процентах: d = 0.298114%
	 РАСЧЕТНЫЕ КОЭФФИЦИЕНТЫ УРАВНЕНИЯ РЕГРЕССИИ:
	 b0 = 6.48866
	 b1 = -0.159952
	 b2 = -0.0359629
	 b3 = -0.376992
	 Остаточная сумма квадратов: Q = 9.89971
	 Число итераций по условию проверки погрешности: n = 184

	 1) Символьное уравнение регрессии:
		 y = b0 + b1*x1 + b2*x2 + b3*x3
	 2) Уравнение регрессии с числовыми оценочными параметрами:
y = 6.48866 + (-0.159952)*x1 + (-0.0359629)*x2 + (-0.376992)*x3
--------------------------------------------------------------

Задание 1

В программе объясните цикл проверки условия погрешности по выходу.

Проверьте программу для случая задания абсолютной погрешности по выходу.

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

Расчет вектора оценочных параметров произведите по встроенной функции $$regress$$ системы MATLAB (см. $$help\mbox{ }regress$$ ).

Расчет вектора оценочных параметров произведите по встроенной функции $$robustfit$$ системы MATLAB (см. $$help\mbox{ }robustfit$$ ). Для функции $$robustfit$$ нет необходимости формировать столбец фиктивной переменной, равной единице.

Расчет вектора оценочных параметров произведите по встроенной функции $$glmfit$$ системы MATLAB (см. $$help\mbox{ }glmfit$$ ). Произвести также расчет значений функции отклика с помощью функции $$glmval$$ системы MATLAB (см. $$help\mbox{ }glmval$$ ).

Примечание. Описание функций $$robustfit$$, $$glmfit$$, $$glmval$$ на русском языке можно найти на сайте http://matlab.exponenta.ru в разделе Statistics Toolbox.

clear,clc
try
   global hep
   delete(hep);
end
 
options.Resize = 'on';
options.WindowStyle = 'normal';
options.Interpreter = 'tex';
f = inputdlg({sprintf('\\bfВектор базисных функций:'),...
sprintf('\\bfОбласть определения:'),...
sprintf('\\bfЧисловые значения параметров b0,b1,...'),...
sprintf('\\bfЗначения аргумента (не менее 4 значений):...')},...
'ИСХОДНЫЕ ДАННЫЕ',1,...
{'x;x^3;sqrt(x)','[0 2]','-5;1.2;3.2;-12.4','0,0.12,0.7,1.34,0.45,1.678 '},options)

if isempty(f)
errordlg('Возможно, нажата клавиша ''Cansel'' ',...
            'Ошибка');
        return    
end
y1 = sym(char(f(1)));
x = str2num(char(f(2)));
b = str2num(char(f(3)));
es = str2num(char(f(4)));
  Y2 = char(y1);
q = find(Y2 == 'x');
if length(q) + 1 ~= length(str2num(char(f(3))) )
        errordlg('Число введенных параметров не соответствует числу параметров введенной функции',...
            'Ошибка ввода параметров');
        return
    end
%---------------- Общее число точек ---------------------
    if length(es) == 1
        if es < length(b)
             errordlg(sprintf('Общее число точек должно быть не менее числа %d - числа параметров',length(b)),...
                 'Ошибка ввода числа точек функции');
             return
         end
          fprintf('\n ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЗАДАННОЙ ФУНКЦИИ\n');
 fprintf('\t Вид параметрических функций:\n\t %s\n',Y2);
fprintf('\t Область определения функции: [%g %g]\n',x(1),x(2));
fprintf('\t Значения введенных параметров:\n');

for J = 1:length(b)
fprintf('\t%g\t',b(J)');
end
fprintf('\n');
disp('-----------------------------------------------------')
er1 = inline(sprintf('[1,%s]',char(f(1))));
erV1 = vectorize(inline(sprintf('[1,%s]',char(f(1)))));
erV = vectorize(inline(sprintf('[%s]',char(f(1)))));
con = char(erV);
n = find(con == ';');
con(n) = ',';
con1 = char(er1);
n1 = find(con1 == ';');
con1(n1) = ',';
symV = sym(con1);
wone = inline(con);
V =x(1) + (x(2) - x(1))*rand(es,1);
fprintf('\t Значения аргумента и параметрических функций:\n');
 
X = [ones(length(V),1),wone(V)];
disp([V,X])
XX = [];
for I = 1:length(V)
    XX = [XX;b'.*X(I,:)];
end
Y = sum(XX,2);
fprintf('\t Значения аргумента и заданной функции:\n');
disp([V,Y])
B = regress(Y,X);
fprintf('\t Расчетные значения параметров заданной функции:\n');
for J = 1:length(B)
    fprintf('\t b%d = %g\n',J - 1, B(J));
end
y = [];
for J = 1:length(symV)
    if J < length(symV)
        if J == 1
         y = [y,sprintf('b%d%s + ',J - 1,'')];
        else
    y = [y,sprintf('b%d*%s + ',J - 1,char(symV(J)))];
end
elseif J == length(symV)
    y = [y,sprintf('b%d*%s',J - 1,char(symV(end)))];
end
end
fprintf('\n\t Заданный вид функциональной зависимости:\n\t y = %s\n',y)
helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'),' ')
%----------------- Несколько точек ------------------------------
elseif length(es) >= length(b)
    if min(es) < x(1) | max(es) > x(2)
    errordlg(sprintf('Значения аргумента функции выходят за область определения [%g, %g]',x(1),x(2)),...
            'Ошибка ввода значений аргумента функции');
        return
    end
     fprintf('\n ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЗАДАННОЙ ФУНКЦИИ\n');
 fprintf(' Вид параметрических функций:\n\t %s\n',Y2);
fprintf(' Область определения функции: [%g %g]\n',x(1),x(2));
disp('-------------------------------------------------------------')
     er1 = inline(sprintf('[1,%s]',char(f(1))));
erV1 = vectorize(inline(sprintf('[1,%s]',char(f(1)))));
erV = vectorize(inline(sprintf('[%s]',char(f(1)))));
con = char(erV);
n = find(con == ';');
con(n) = ',';
con1 = char(er1);
n1 = find(con1 == ';');
con1(n1) = ',';
symV = sym(con1);
wone = inline(con);
V = es';
fprintf(' Значения аргумента и значения параметрических функций:\n');
X = [ones(length(V),1),wone(V)];
disp([V,X])
XX = [];
for I = 1:length(V)
    XX = [XX;b'.*X(I,:)];
end
Y = sum(XX,2);
fprintf(' Значения аргумента и значения заданной функции:\n');
disp([V,Y])
B = regress(Y,X);
fprintf(' Расчетные значения параметров заданной функции:\n');
for J = 1:length(B)
    fprintf('\t b%d = %g\n',J - 1,B(J));
end

y = [];
for J = 1:length(symV)
    if J < length(symV)
        if J == 1
         y = [y,sprintf('b%d%s + ',J - 1,'')];

        else
    y = [y,sprintf('b%d*%s + ',J - 1,char(symV(J)))];
end

elseif J == length(symV)
    y = [y,sprintf('b%d*%s',J - 1,char(symV(end)))];
end
end
fprintf('\n Заданный вид функциональной зависимости:\n\t y = %s\n',y);

global hep
hep = helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'), '');
else 
    errordlg(sprintf('Необходимо общее число аргументов не менее %d',length(b)), 'Ошибка!')
    return
end

clear,clc,close all
x = ... [0;5;10;15;20;25;30;35;40;45;50;55;60;65;70;72;74;75;76;78;80;82;84;85;86;88;90];
f = ...
[149.0;146.3;145.8;149.8;154.9;164.7;178.2;196.2;221.8;259.8;299.2;331.5;338.1;304.4;238.5;205.3;170.9;153.9;
  136.6;103.9;74.2;48.6;29.7;19.9;12.4;3.1;0];   
n = length(x); %% length(x) == length(f)
 Lmin = f'*f; %% Квадратичный критерий - сумма квадратов
line(x,f,'linew', 2)
grid on
xlabel('\bf - - - - - - - - - x - - - - - - - - -');
ylabel('\bf f(x)');
 qprob = 0.4; % Уровень значимости
%%%% Выбор степени полинома
for m = 1 : n
    psi(:, m) = x.^(m-1); % psi - Специальная функция (полигамма)
    if m > 1
       for k = 1 : m-1
          psi(:, m) = psi(:, m) - ((psi(:,m))'*psi(:,k))*psi(:, k);
       end
    end
    psi(:, m) = psi(:,m)/norm(psi(:,m)); % нормировка
 
    bk = f'*psi(:,m);
    DOLD = Lmin/(n-m); %%%Старая дисперсия
    Lmin = Lmin - bk^2;
    Dnew = Lmin/(n-m-1); %%% Новая дисперсия
% fprintf('\n');
    if m > 1 % проверка степени начиная с 2
       if Dnew/DOLD < finv(qprob, n-m, n-m-1) 
% fprintf(' Степень %d следует учитывать\n', m-1);
           stp = m-1;
       else
% fprintf(' Степень %d не нужно учитывать\n', m-1);
 
if m > 2
 stp = m-2;
break;
else
  stp = m - 1;  
break;
end
        end
    end
end
 PO = polyfit(x, f, stp); 
fprintf('\n Коэффициенты аппроксимируещего полинома:\n'); 
for J = 1 : length(PO)
fprintf('%15g\n', PO(J));
end

fprintf('\n');
 yPO = polyval(PO, x);
 
line(x, yPO, 'marker', 'o', 'color','r')
title(sprintf('%s Аппроксимация экспериментальной кривой полиномом %d степени','\bf', stp));
legend('\bfЭксперимент','\bfАппроксимация', 'location', 'best')

clear,clc
%---------------------------------------------------------------
options.Resize = 'on';
   options.WindowStyle = 'normal';
   options.Interpreter = 'tex';
   yx = inputdlg({'\bfСредние значения функций отклика:',...
        '\bfСредние значения аргументов функции отклика:',...
        '\bfЧисло измерений:',...
        '\bfВеличина стандартного отклонения функции откли-ка:.....'},...
'ИСХОДНЫЕ ДАННЫЕ',1,{'1,2,3,4,5,6,7,8','11,23,4,7,9','50','0.1'},options);

% БАЗОВЫЕ ВЕКТОРЫ ВЫХОДА, ВХОДА, ЧИСЛО ИЗМЕРЕНИЙ, 
% СТАНДАРТНОЕ ОТКЛОНЕНИЕ
Yb = str2num(char(yx(1)));
Ub = str2num(char(yx(2)));
r = str2num(char(yx(3)));
s = str2num(char(yx(4)));
%-------------------------------------------------------------------
% ПРОВЕРКА ВВОДА ДАННЫХ
if ( isempty(Yb) == 1 | isreal(Yb) == 0 | isinf(Yb) == 1 )
    errordlg(sprintf('Значения функций должны составлять вектор действительных чисел'),'Ошибка ввода значений функций')
   return

elseif ( min(size(Yb)) > 1  max(size(Yb)) > 1 )
   errordlg(sprintf('Значения функций должны составлять вектор, не матрицу'),'Ошибка ввода значений функций')
   return
end
%----------------------------------------------------------------
%% ПРИВЕДЕНИЕ К ВЕКТОРУ-СТРОКЕ
[ry,vy] = size(Yb);
if ( ry == 1 )
    Yb = Yb;
elseif ( ry > 1 )
    Yb = Yb';
end
if ( isempty(Ub) | ~isreal(Ub) | isinf(Ub) )
    errordlg(sprintf('Значения аргументов должны составлять вектор действительных чисел'), ...
        'Ошибка ввода значений аргументов')
    return
elseif ( min(size(Ub)) > 1  max(size(Ub)) > 1 )
   errordlg('Значения аргументов функций должны составлять вектор, не матрицу','Ошибка ввода значений аргументов функций')
   return
elseif ( isempty(find(Ub == 0)) == 0 )
    errordlg(sprintf('Значения аргументов функций не должны быть нулевыми'),'Ошибка ввода значений аргументов')
   return
end
%--------------------------------------------------------------
%% ПРИВЕДЕНИЕ К ВЕКТОРУ-СТРОКЕ
[ru,vu] = size(Ub);
if ( ru == 1 )
    Ub = Ub;
elseif ( ru > 1 )
    Ub = Ub';
end
if ( isempty(r) | ~isreal(r) | isinf(r) | r <= 0 )
    errordlg('Число измерений должно быть натуральным числом',...
               'Ошибка ввода числа измерений')
    return
elseif ( length(r) > 1 | r ~= round(r))
   errordlg('Число измерений должно быть одним натуральным числом',...
             'Ошибка ввода числа измерений')
   return

elseif ( r <= length(Yb)*(length(Ub) + 1) )
    errordlg(sprintf('Число измерений должно быть не меньше произведения числа функций на число аргументов плюс 1: больше %d',...
length(Yb)*(length(Ub) + 1)), 'Ошибка ввода числа измерений')
    return
end
%-----------------------------------------------------------------
if ( isempty(s) | ~isreal(s) | isinf(s) )
errordlg('Стандартное отклонение должно быть действительным числом',...
 'Ошибка ввода стандартного отклонения')
    return
elseif ( length(s) > 1 )
   errordlg('Стандартное отклонение должно быть одним действительным числом','Ошибка ввода стандартного отклонения')
   return
elseif ( s <= 0 )
    errordlg('Стандартное отклонение должно быть больше нуля',...
        'Ошибка ввода отклонения')
    return
end
%----------------------------------------------------------------
% % ЗАДАНИЕ r  ИЗМЕРЕНИЙ ФУНКЦИЙ ОТКЛИКА
Y = [];
 for k = 1 : length(Yb)
    for I = 1 : r
         Y(:,k) = Yb(k) + normrnd(0,s,r,1);
    end
end
% % ЗАДАНИЕ r  ИЗМЕРЕНИЙ АРГУМЕНТОВ ФУНКЦИИ ОТКЛИКА 
 X = [];
for m = 1 : length(Ub)
    for J = 1 : r
    X(:,m) = Ub(m) + rand(r,1);
end
end
%% ВВЕДЕНИЕ ФИКТИВНОГО АРГУМЕНТА, РАВНОГО ЕДИНИЦЕ
X1 = [ones(r,1),X];
%% РАСЧЕТ ПАРАМЕТРОВ ФУНКЦИИ ОТКЛИКА
 B = (inv(X1'*X1)*X1'*Y)';
 %% ВЫВОД РЕЗУЛЬТАТОВ
 fprintf('\n\t ИСХОДНЫЕ ДАННЫЕ ЭКСПЕРИМЕНТА:\n')
fprintf('\t Вектор-строка значений функции отклика:\n')
fprintf('\t %g',Yb)
fprintf('\n')
fprintf('\t Вектор-строка значений аргументов:\n')
fprintf('\t %g',Ub)
fprintf('\n')
fprintf('\t Заданное число повторных измерений: r = %d\n',r)
fprintf('\t Стандартное отклонение повторных измерений функций: s = %g\n',s)
disp('-----------------------------------------------------------')
fprintf('\t РЕЗУЛЬТАТЫ ПАССИВНОГО МНОГОМЕРНОГО ЭКСПЕРИМЕНТА\n')
fprintf('\t ОЦЕНКА МАТРИЦЫ ПАРАМЕТРОВ:\n')
 disp(B)
 %% РАСЧЕТНЫЕ ЗНАЧЕНИЯ ФУНКЦИЙ
 
Ybpac = B*[1,Ub]';
disp('------------------------------------------------------------')
fprintf('\t СРАВНЕНИЕ ЗНАЧЕНИЙ ЗАДАННЫХ ФУНКЦИЙ И РАСЧЕТНЫХ:\n')
disp([Yb' Ybpac])
disp('------------------------------------------------------------')
helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'),'')

ИСХОДНЫЕ ДАННЫЕ ЭКСПЕРИМЕНТА
	 Вектор-строка значений функции отклика:
	 1	 2	 3	 4	 5	 6	 7	 8
	 Вектор-строка значений аргументов:
	 11	 23	 4	 7	 9
	 Заданное число повторных измерений: r = 50
	 Стандартное отклонение повторных измерений функций: s = 0.1
-----------------------------------------------------------
	 РЕЗУЛЬТАТЫ ПАССИВНОГО МНОГОМЕРНОГО ЭКСПЕРИМЕНТА
	 ОЦЕНКА МАТРИЦЫ ПАРАМЕТРОВ:
   -1.5280   -0.0296    0.0969   -0.0103   -0.0108    0.0767
    3.1418   -0.0457   -0.0015    0.0794   -0.0597   -0.0519
    2.2212   -0.0010    0.0370   -0.0022    0.0191   -0.0225
    4.3340   -0.1121    0.0989   -0.0480    0.0087   -0.1287
    1.3881    0.0194    0.0830   -0.0053    0.0519    0.1112
    8.1169   -0.0406   -0.1085    0.0382   -0.0001    0.0760
    9.9112   -0.0201   -0.0095    0.0651   -0.0843   -0.0143
    8.1976    0.0707   -0.0004   -0.0371   -0.0408   -0.0563
------------------------------------------------------------
	 СРАВНЕНИЕ ЗНАЧЕНИЙ ЗАДАННЫХ ФУНКЦИЙ И РАСЧЕТНЫХ:
    1.0000    0.9483
    2.0000    2.0366
    3.0000    2.9832
    4.0000    4.0865
    5.0000    4.8532
    6.0000    6.0105
    9.0000    9.0124
    8.0000    8.0252
------------------------------------------------------------