mirror of
https://github.com/eddyem/lectures.git
synced 2025-12-06 10:45:09 +03:00
add 03-pract
This commit is contained in:
Edward Emelianov
parent
0d8fbff1dd
commit
20478a2c90
BIN
Komp_obr_SFedU/03_Pract.pdf
Normal file
BIN
Komp_obr_SFedU/03_Pract.pdf
Normal file
Binary file not shown.
569
Komp_obr_SFedU/03_Pract.tex
Normal file
569
Komp_obr_SFedU/03_Pract.tex
Normal file
@ -0,0 +1,569 @@
|
||||
\documentclass[a4paper,12pt]{extarticle}
|
||||
\usepackage{/home/eddy/ed, verbatim}
|
||||
\title{Практикум \No3: погрешности, метод наименьших квадратов}
|
||||
\author{}\date{}\nocolon
|
||||
|
||||
\long\def\task#1{\noindent\leavevmode\refstepcounter{sect}\llap{\textbf{\thesect}\;}\indent\textit{#1}}
|
||||
\def\t#1{{\upshape\ttfamily #1}}
|
||||
|
||||
\begin{document}
|
||||
\maketitle
|
||||
\section{Погрешности}
|
||||
\subsection{}
|
||||
Найдем общую среднюю совокупности, состоящей из следующих трех групп:
|
||||
$$\renewcommand{\arraystretch}{0}
|
||||
\begin{tabular}{|r|c|c|c|c|c|c|}
|
||||
\hbox to 5cm{}&\hbox to 1.3cm{}& \hbox to 1.3cm{}& \hbox to 1.3cm{}&
|
||||
\hbox to 1.3cm{}& \hbox to 1.3cm{}& \hbox to 1.3cm{}\\
|
||||
\hline
|
||||
\strut Группа& \multicolumn{2}{|c|}{I} & \multicolumn{2}{|c|}{II} &
|
||||
\multicolumn{2}{|c|}{III} \\
|
||||
\hline
|
||||
\strut Значение признака&1&3&2&4
|
||||
&3&6\\
|
||||
\strut Частота признака&11&34&22&28&31&14\\\hline
|
||||
\strut Объем выборки&\multicolumn{2}{|c|}{$11+34=45$}&\multicolumn{2}{|c|}{$22+28=50$}&
|
||||
\multicolumn{2}{|c|}{$31+14=45$}\\
|
||||
\hline
|
||||
\end{tabular}
|
||||
$$
|
||||
|
||||
Для начала найдем групповые средние: $\aver{x_1}$, $\aver{x_2}$ и~$\aver{x_3}$:
|
||||
\begin{verbatim}
|
||||
x1 = (11*1 + 34*3)/45
|
||||
x2 = (22*2 + 28*4)/50
|
||||
x3 = (31*3 + 14*6)/45
|
||||
\end{verbatim}
|
||||
|
||||
Теперь найдем их среднее:
|
||||
\begin{verbatim}
|
||||
X = (x1*45 + x2*50 + x3*45)/(45 + 50 + 45)
|
||||
\end{verbatim}
|
||||
|
||||
Однако, при работе с большими массивами данных лучше использовать преимущества
|
||||
матричной алгебры:
|
||||
\begin{verbatim}
|
||||
xi = [1 3 2 4 3 6];
|
||||
ni = [11 34 22 28 31 14 ];
|
||||
N = sum(ni)
|
||||
X = sum(xi.*ni/N)
|
||||
\end{verbatim}
|
||||
|
||||
Найдем\к генеральную дисперсию\н и генеральное среднеквадратичное отклонение
|
||||
данной выборки:
|
||||
\begin{verbatim}
|
||||
D = sum(ni.*(xi-X).^2)/N
|
||||
sigma=sqrt(D)
|
||||
\end{verbatim}
|
||||
Кроме того, определить среднеквадратичное отклонение ряда~$x$ можно при помощи
|
||||
команды \verb'std(x)'.
|
||||
|
||||
\subsection{}
|
||||
Рассмотрим ряд измерений некоторой физической величины~$x$. Результаты
|
||||
серии измерений заданы таблицей ($\nu_i$~-- частота соответствующего значения~$x_i$):
|
||||
$$
|
||||
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}
|
||||
\hline
|
||||
$x_i$& 31 & 28 & 34 & 26 & 35 & 30 & 34 & 32 & 40 & 20 \\
|
||||
$\nu_i$& 20 & 12 & 10 & 5 & 7 & 20 & 12 & 19 & 4 & 2 \\
|
||||
\hline
|
||||
\end{tabular}
|
||||
$$
|
||||
Известно, что некоторые результаты могут быть заведомо ошибочными.
|
||||
Нам необходимо оценить среднее значение данной величины, исключив ошибочные
|
||||
результаты.
|
||||
Составим массивы величины~$x$ и соответствующих частот~$n$:
|
||||
\begin{verbatim}
|
||||
x = [ 31 28 34 26 35 30 34 32 40 20];
|
||||
n = [ 20 12 10 5 7 20 12 19 4 2];
|
||||
\end{verbatim}
|
||||
Отобразив данные на графике (\verb"plot(x,n,'o')") можно заметить, что действительно
|
||||
некоторые значения сильно отклоняются от положения, которое они занимали бы при
|
||||
нормальном распределении.
|
||||
|
||||
Найдем среднее значение величины~$x$ и ее среднеквадратичное отклонение:
|
||||
\begin{verbatim}
|
||||
X = sum(x.*n)/sum(n)
|
||||
sigma = sqrt(sum(n.*(x-X).^2)/sum(n))
|
||||
\end{verbatim}
|
||||
Определим границы доверительного интервала $[a,b]$ в пределах трех~$\sigma$:
|
||||
\begin{verbatim}
|
||||
a = X-3*sigma
|
||||
b = X+3*sigma
|
||||
\end{verbatim}
|
||||
Теперь исключим из выборки значения, выходящие за пределы интервала. При помощи
|
||||
функции \verb'find' можно найти индексы членов массива, удовлетворяющих заданному
|
||||
условию. Исключить лишние элементы можно так:
|
||||
\begin{verbatim}
|
||||
idx = find(x < a);
|
||||
x(idx) = []; n(idx) = [];
|
||||
idx = find(x > b);
|
||||
x(idx) = []; n(idx) = [];
|
||||
\end{verbatim}
|
||||
|
||||
Теперь повторим вычисление \verb'X' и \verb'sigma', \verb'a' и \verb'b'.\к
|
||||
Для того, чтобы вызвать из истории команд строку, начинающуюся с определенных
|
||||
символов, наберите один--два первых символа и нажмите клавишу <<вверх>>\н.
|
||||
Таким образом можно быстро вызвать из истории команд нужную вам команду, не
|
||||
перебирая все промежуточные.
|
||||
|
||||
Теперь проверим, не влияет ли <<подозрительное>> значение $x=40$ на точность
|
||||
измерения. Найдем медиану нашего ряда и оценим доверительный интервал по медиане.
|
||||
Для этого нам необходимо построить новый вектор~\verb'newx', в котором значения
|
||||
величины~$x$ будут содержаться столько раз, какова их частота:
|
||||
\begin{verbatim}
|
||||
newx = []; for a = [1:length(n)]
|
||||
newx = [newx ones(1,n(a)).*x(a)];
|
||||
endfor
|
||||
med = median(newx)
|
||||
a = med-3*sigma
|
||||
b = med+3*sigma
|
||||
\end{verbatim}
|
||||
Действительно, значение $x=40$ выбивается из доверительного интервала.
|
||||
Удалим его:
|
||||
\begin{verbatim}
|
||||
idx = find(x==40);
|
||||
x(idx) = []; n(idx) = [];
|
||||
\end{verbatim}
|
||||
и найдем~$\aver{x}$, близкое к истинному:
|
||||
\begin{verbatim}
|
||||
X = sum(x.*n)/sum(n)
|
||||
sigma = sqrt(sum(n.*(x-X).^2)/sum(n))
|
||||
a = X-3*sigma
|
||||
b = X+3*sigma
|
||||
find(x>b)
|
||||
find(x<a)
|
||||
\end{verbatim}
|
||||
Итак, все оставшиеся значения~$x_i$ удовлетворяют критерию
|
||||
<<трех сигм>>, следовательно, можно записать ответ: $x=31.3\pm2.3$.
|
||||
|
||||
\subsection{}
|
||||
Теперь определим доверительный интервал величины $\aver{x}$ с
|
||||
надежностью~95\% при помощи распределения Стьюдента. Для этого в Octave существует
|
||||
функция~\verb'ttest'. В простейшем случае вида \verb'h=ttest(x)' она возвращает
|
||||
вероятность отклонения гипотезы о нормальном распределении величины~$x$ с
|
||||
математическим ожиданием $\mean{x}=0$. Проверка даст результат:~1. Действительно,
|
||||
математическое ожидание нашей величины далеко не равно нулю. Второй аргумент
|
||||
функции \verb'ttest' задает предполагаемое математическое ожидание. Проверим:
|
||||
\verb'h=ttest(x,X)'. Получаем: \verb'h=0'. Т.е., можно принять гипотезу
|
||||
о гауссовой форме распределения величины~$x$ около ее среднего значения.
|
||||
Оценить 95\%-й доверительный интервал величины~$x$ можно при помощи расширенного
|
||||
вывода функции \verb'ttest' в форме \verb'[h,p,ci]=ttest(x,X)'. В этом случае
|
||||
параметр \verb'h' сообщает о степени ненадежности гипотезы, \verb'p' равен
|
||||
вероятности совпадения величины~\verb'X' с математическим ожиданием ряда~\verb'x',
|
||||
\verb'ci' сообщает границы 95\%-го доверительного интервала.
|
||||
Определим доверительный интервал для нашего ряда без исключения заведомо ложных
|
||||
результатов и с их исключением:
|
||||
\begin{verbatim}
|
||||
x = [ 31 28 34 26 35 30 34 32 40 20];
|
||||
n = [ 20 12 10 5 7 20 12 19 4 2];
|
||||
newx = []; for a = [1:length(n)]
|
||||
newx = [newx ones(1,n(a)).*x(a)];
|
||||
endfor
|
||||
|
||||
[h,p,ci] = ttest(newx, X)
|
||||
newx = []; for a = 1:8
|
||||
newx = [newx ones(1,n(a)).*x(a)];
|
||||
endfor
|
||||
[h,p,ci] = ttest(newx, X)
|
||||
\end{verbatim}
|
||||
|
||||
Итак, в обоих случаях гипотеза о соответствии распределения величины~$x$
|
||||
нормальному распределению принимается, однако, во втором случае вероятность
|
||||
определения математического ожидания~$\mean{x}$ выше, и доверительный интервал
|
||||
уже, что явно свидетельствует о большей надежности вычислений.
|
||||
|
||||
\subsection{}
|
||||
Octave предоставляет огромный набор инструментальных средств. Однако, при работе с большим
|
||||
количеством однообразных данных приходится много раз повторять одни и те же команды. Эту задачу
|
||||
можно упростить, создав\ж скрипт\н (или m-файл). Скрипт представляет собой описание и реализацию
|
||||
пользовательской функции, которая вызывается из командной строки Octave аналогично
|
||||
любой команде, однако может содержать значительное количество инструкций,
|
||||
облегчающих работу пользователя.
|
||||
|
||||
M-файл может содержать любые инструкции. Если он не начинается со слова
|
||||
\verb'function', выполняется все его содержимое. Удобнее,
|
||||
однако, создать m-файл в виде функции, принимающей в качестве аргументов
|
||||
необходимые переменные и возвращающей определенные величины.
|
||||
|
||||
Заголовок файла функции имеет вид
|
||||
\begin{verbatim}
|
||||
%
|
||||
% Комментарий, отображающийся при введении команды help имя_функции
|
||||
%
|
||||
function [возвращаемые величины] = имя_функции(входные, аргументы)
|
||||
\end{verbatim}
|
||||
|
||||
Далее следуют операторы, выполняемые в теле функции. Если после команды вы
|
||||
пропустите символ точки с запятой, ее вывод будет отображен на экране.
|
||||
|
||||
Итак, создадим m-файл, осуществляющий проверку выборки на корректность
|
||||
при помощи критерия <<трех сигм>>:
|
||||
\verb'three_s.m'.
|
||||
{\small
|
||||
\verbatiminput{Materials4Pract/03/three_s.m}
|
||||
}
|
||||
|
||||
При запуске скрипта по умолчанию Octave его ищет в текущей директории. Однако, можно добавить любую
|
||||
директорию со скриптами в список поиска (\t{path}) при помощи команды \t{addpath}.
|
||||
|
||||
Запустить данный скрипт можно командой \verb'[X sigma] = three_s(x,n)'.
|
||||
|
||||
\subsection{}
|
||||
Зачастую физику-экспериментатору приходится проверять нулевую гипотезу о
|
||||
равенстве средних двух независимых наборов данных.
|
||||
Пусть в результате одного измерения некоторой физической величины~$x$ был
|
||||
получен ряд данных:
|
||||
\begin{verbatim}
|
||||
x1 = [ 47.78 36.40 35.66 8.93 40.42 54.16 51.76 44.32 46.19 50.75];
|
||||
\end{verbatim}
|
||||
Затем было произведено независимое измерение этой же физической величины
|
||||
при других условиях эксперимента. При этом был получен ряд:
|
||||
\begin{verbatim}
|
||||
x2 = [ 44.09 46.75 44.20 7.99 47.74 75.07 62.48 44.43 34.73 55.26];
|
||||
\end{verbatim}
|
||||
Требуется проверить нулевую гипотезу о равенстве математических ожиданий
|
||||
данных величин.
|
||||
|
||||
Для проверки данной гипотезы существует функция Octave \verb'ttest2',
|
||||
\verb'ttest2(x1,x2)' даст ответ: \verb'ans=0', т.е.
|
||||
гипотеза о неравенстве математических ожиданий наших двух рядов отклонена
|
||||
на 95\%-м уровне. Для определения доверительного интервала и вероятности
|
||||
равенства математических ожиданий воспользуемся расширенным выводом команды:
|
||||
\begin{verbatim}
|
||||
[h p ci] = ttest2(x1, x2)
|
||||
\end{verbatim}
|
||||
Таким образом, вероятность того, что математические ожидания выборок равны,
|
||||
составляет лишь~$p=51\%$, при этом доверительный интервал математического
|
||||
ожидания разности $x_1-x_2$ достаточно широк: $c_i=[-19.2,9.9]$, т.е.
|
||||
математические ожидания данных рядов могут разниться на~$4.6$ со среднеквадратичным
|
||||
отклонением~$\sigma=14.6$.
|
||||
|
||||
Большая ширина доверительного интервала говорит о том, что данные в рядах~$x_1$
|
||||
и~$x_2$ получены с низкой надежностью. Однако, найдя медианы рядов~$x_1$,
|
||||
$x_2$ и совмещенного ряда $(x_1;x_2)$ можно попытаться с достаточно высокой
|
||||
степенью вероятности оценить математическое ожидание величины~$x$.
|
||||
|
||||
|
||||
\section{Метод наименьших квадратов}
|
||||
\subsection{}
|
||||
Случайная погрешность физического измерения имеет природу, аналогичную белому шуму,
|
||||
поэтому для начала рассмотрим простейшие методы очистки одномерных сигналов вида
|
||||
$y=y(t)$ от шумов.
|
||||
|
||||
Используем массив из десяти сигналов, зашумленных с одинаковым уровнем~SNR:
|
||||
\begin{verbatim}
|
||||
x=[0:0.05:20];
|
||||
y=sin(x*10).*(0.5+sawtooth(x*pi/5)/2);
|
||||
for a=[1:10]
|
||||
y1(a,:)=awgn(y,1,'measured');
|
||||
endfor
|
||||
\end{verbatim}
|
||||
|
||||
Вид цикла \verb'for' отличается от языков программирования вроде C: цикл
|
||||
поочередно перебирает все значения переменной \verb'a'. Если бы мы заранее
|
||||
инициализировали ее массивом, можно было бы просто написать \verb'for a'.
|
||||
Цикл \verb'for' заканчивается командой \verb'endfor'. Двоеточие в адресации \verb'y(a,:)' означает,
|
||||
что мы выбираем\ж все\н элементы по второй координате (т.е. приравнивание производится к целой
|
||||
строке). Еще одним отличием от языков программирования является динамическое
|
||||
расширение матриц: нет необходимости в начале работы с ней сообщать ее предельный
|
||||
размер.
|
||||
|
||||
Итак, мы получили массив \verb'y1', в строках которого содержатся зашумленные
|
||||
варианты одного и того же сигнала. Можно отобразить их все графически командой
|
||||
\verb'plot(x,y1)', а можно и с оригиналом: \verb'plot(x,y,"linewidth",2, x, y1)'.
|
||||
Оценить зашумленность сигнала можно командой \verb"plot(y,y1,'.')". Если бы
|
||||
сигналы в \verb'y1' совпадали с \verb'y', мы увидели бы отрезок с коэффициентом
|
||||
наклона~1. Чем дальше форма полученной фигуры от такого отрезка, тем больше
|
||||
зашумленность сигнала.
|
||||
|
||||
Для восстановления сигнала из десяти измерений попробуем усреднить наборы сигналов
|
||||
и найти их медиану:
|
||||
\begin{verbatim}
|
||||
y_mean = mean(y1);
|
||||
y_med = median(y1);
|
||||
plot(x,[y;y_mean;y_med]);
|
||||
legend("original", "mean", "median");
|
||||
\end{verbatim}
|
||||
|
||||
Оба восстановленных сигнала имеют примерно одинаковые величины
|
||||
и довольно близки к реальной функции (особенно на участках с большой амплитудой
|
||||
сигнала). Однако, как мы увидим впоследствии, если к сигналу добавлен шум
|
||||
типа <<соль/перец>>, медианная фильтрация будет работать намного эффективнее
|
||||
фильтрации по среднему арифметическому.
|
||||
|
||||
\subsection{}
|
||||
Рассмотрим линейную зависимость $y=ax+b$, заданную таблично в виде~$y=y(x)$.
|
||||
Для определения методом наименьших квадратов коэффициентов линейной (а также высших степеней)
|
||||
зависимости служит функция \verb'polyfit(x,y,n)'. Она содержит три аргумента: $x$~-- вектор
|
||||
аргумента, $y$~-- вектор функции, $n$~-- степень аппроксимирующего полинома. Ее результат в
|
||||
простейшем случае представляет собой вектор коэффициентов (начиная со старшей степени).
|
||||
Если функцию вызвать как \verb'[p,S] = polyfit(x,y,n)', вектор~$p$ будет содержать коэффициенты,
|
||||
а в структуре~$S$ будут содержаться такие данные, как степени свободы~(df) и норма отклонений
|
||||
данных от аппроксимирующей кривой (normr). Для восстановления полученной зависимости
|
||||
используется функция \verb'polyval(p,x)', где $p$~-- полученный функцией \verb'polyfit'
|
||||
вектор коэффициентов, $x$~-- вектор аргумента. В таком виде функция возвращает вектор
|
||||
восстановленной функции. В виде \verb'[y, delta] = polyval(p,x,S)' функция возвращает массив
|
||||
погрешностей (т.е. в каждой точке восстановленные значения функции можно представить в
|
||||
виде $y=y\pm delta$, т.е. оценить абсолютную погрешность восстановления можно при помощи
|
||||
команды \verb'mean(delta)'.
|
||||
|
||||
Найдем коэффициенты модельной зависимости. Пусть $y=7.15x+4.22$. Построим
|
||||
векторы, соответствующие аргументу и функции:
|
||||
$$
|
||||
\verb'x = [0:100]; y = 7.15*x + 4.22;'
|
||||
$$
|
||||
Зашумим сигнал для получения разброса точек~$y_i$:
|
||||
$$
|
||||
\verb"y1 = awgn(y,10,'measured');"
|
||||
$$
|
||||
Отобразим на экране оба ряда: \verb"plot(x,y,x,y1,'o')" (запись \verb"'o'"
|
||||
означает, что график будет отображаться кружк\'ами). Разброс данных
|
||||
достаточно велик. Определим коэффициент корреляции: \verb'corr(x,y1)'.
|
||||
Он довольно близок к единице, следовательно, мы можем попытаться получить
|
||||
коэффициенты линейной зависимости и восстановить функцию:
|
||||
\begin{verbatim}
|
||||
[p,S] = polyfit(x,y1,1); % коэффициенты a и b
|
||||
[y2, delta] = polyval(p,x,S); % восстановленный вектор
|
||||
plot(x,y1,'o',x,[y;y2]) % все три графика
|
||||
legend("noicy", "original", "fitted");
|
||||
mean(delta) % абсолютная ошибка
|
||||
mean(delta)/mean(y) % относительная ошибка
|
||||
\end{verbatim}
|
||||
|
||||
\subsection{}
|
||||
Можно найти приближение методом наименьших квадратов и другим способом. Пусть
|
||||
$Y$~-- вектор--столбец значений функции, $A=(a,b)^{\rm Tr}$~-- вектор--столбец
|
||||
коэффициентов
|
||||
разложения. Тогда условие $y_i=ax_i+b$ можно представить в виде матричного
|
||||
произведения $Y=XA$. Второй столбец матрицы~$X$ целиком состоит из единиц, а
|
||||
в первом находится последовательность значений~$x_i$. В этом случае нахождение
|
||||
коэффициентов сводится к решению системы линейных уравнений $y_i=ax_i+b$,
|
||||
дающему минимальную невязку. Такое решение находится при помощи операции
|
||||
левостороннего матричного деления: $X\backslash Y$. Решим предыдущий пример таким
|
||||
способом.
|
||||
\begin{verbatim}
|
||||
X = [x' ones(size(x'))]; % создаем матрицу аргумента
|
||||
% (т.к. x и y1 - строки, транспонируем их)
|
||||
A = X\y1' % находим коэффициенты
|
||||
% и отображаем их на экране
|
||||
\end{verbatim}
|
||||
|
||||
Полученные значения должны быть примерно равны найденным предыдущим способом.
|
||||
Как мы увидим далее, такой способ нахождения корней аппроксимации пригоден не
|
||||
только для полиномиальных, но и для многих других функций.
|
||||
|
||||
\subsection{}
|
||||
Попробуем создать квадратичную зависимость и аппроксимировать ее методом наименьших квадратов.
|
||||
Пусть зависимость на отрезке $[0,100]$ имеет вид $y=2.4x^2-0.87x+2.13$. Создадим соответствующие
|
||||
массивы данных, добавим шум с SNR=20\,дБ и отобразим оба сигнала на графике:
|
||||
\begin{verbatim}
|
||||
x = [1:100];
|
||||
y = 2.4*x.^2-0.87*x+2.13;
|
||||
y1 = awgn(y,20,'measured');
|
||||
plot(x,[y;y1]);
|
||||
\end{verbatim}
|
||||
|
||||
Теперь создадим вектор коэффициентов аппроксимации полиномом второй степени
|
||||
восстановим функцию и отобразим на графике:
|
||||
\begin{verbatim}
|
||||
[p, S] = polyfit(x, y1, 2);
|
||||
[y2, DELTA] = polyval(p, x, S);
|
||||
plot(x,[y;y2]);
|
||||
legend("original", "restored")
|
||||
\end{verbatim}
|
||||
|
||||
Отобразим на экране найденные коэффициенты:
|
||||
\verb'p'. Рассчитаем среднее квадратичное отклонение аппроксимации (\verb'mean(DELTA)').
|
||||
Также рассчитаем относительную ошибку аппроксимации \verb'mean(DELTA)/mean(y1)'.
|
||||
|
||||
И второй способ:
|
||||
\begin{verbatim}
|
||||
X = [(x.^2)' x' ones(size(x'))];
|
||||
A = X\y1'
|
||||
\end{verbatim}
|
||||
|
||||
\subsection{}
|
||||
Однако, чаще всего функциональные зависимости имеют иные виды зависимости. Допустим, нам
|
||||
известно, что измеряемая величина изменяется по закону
|
||||
\begin{equation}
|
||||
y=a_0+a_1\e^{-t}+a_2te^{-t}.
|
||||
\label{exp_y}
|
||||
\end{equation}
|
||||
Для аппроксимации такой функцией можно представить уравнение~\eqref{exp_y} в матричном
|
||||
виде $Y=TA$, где $T$~-- функциональная матрица, у которой в первом столбце
|
||||
размещены единицы (соответствует умножению на~$a_0$), во втором~--- функция
|
||||
$\e^{-t}$, а в третьем~--- $t\e^{-t}$. Найти коэффициенты~$A$ можно при помощи
|
||||
оператора левого деления: $A=T\backslash Y$.
|
||||
\begin{verbatim}
|
||||
t = [0 0.3 0.8 1.1 1.6 2.3]'; % сразу вводим данные в столбцах
|
||||
y = [0.6 0.67 1.01 1.35 1.47 1.25]';
|
||||
T = [ones(size(t)) exp(-t) t.*exp(-t)];
|
||||
A = T\y
|
||||
\end{verbatim}
|
||||
Теперь отобразим данные на графике:
|
||||
\begin{verbatim}
|
||||
x = [0:0.1:2.5]';
|
||||
Y = [ones(size(x)) exp(-x) x.*exp(-x)]*A;
|
||||
plot(x,Y, t,y,'o')
|
||||
\end{verbatim}
|
||||
|
||||
\subsection{}
|
||||
Для коррекции наведения и сопровождения телескопов используется СКН~--- система коррекции
|
||||
наведения, которая учитывает различного рода ошибки (гнутие осей и неперпендикулярность осей и
|
||||
т.п.). У БТА данные ошибки выражаются полиномами:
|
||||
$$dA = K_0 + K_1\frac{1}{\tg Z} + K_2\frac{1}{\sin Z} - K_3\frac{\sin A}{\tg Z}
|
||||
+K_4\frac{\cos\delta\cos P}{\sin Z},$$
|
||||
$$dZ = K_5 + K_6\sin Z + K_7\cos Z + K_3\cos A + K_4 \cos\phi\sin A.$$
|
||||
|
||||
Здесь:
|
||||
\begin{description}
|
||||
\item[$dA$, $dZ$] погрешности наведения по азимуту и зенитному расстоянию;
|
||||
\item[$K_x$] коэффициенты ошибок;
|
||||
\item[$\phi$] широта места наблюдения;
|
||||
\item[$t$] часовой угол;
|
||||
\item[$P$] параллактический угол.
|
||||
\end{description}
|
||||
|
||||
Для получения этих коэффициентов необходимо провести наблюдение в нескольких десятках равномерно
|
||||
распределенных по небесной полусфере точек (за исключением областей запрета, $Z<5\degr$, и около
|
||||
горизонта, $Z>70\degr$). Далее в каждом поле вычисляется астрометрия и определяется погрешность
|
||||
наведения. Составляется таблица (например, \t{2015\_09\_30\_pf.tab}, по которой и необходимо
|
||||
вычислить коэффициенты. Вычислять коэффициенты будем следующим скриптом:
|
||||
|
||||
\verbatiminput{Materials4Pract/03/getSKNcoeff.m}
|
||||
|
||||
Запускаем: \t{SKN = getSKNcoeff('2015\_09\_30\_pf.tab')}. Строятся графики остаточных невязок и
|
||||
выводятся значения всех коэффициентов.
|
||||
|
||||
Из-за линейной зависимости коэффициентов задача их вычисления является некорректной, а ввиду
|
||||
малости объема экспериментального материала, решать задачу будем итерациями, на каждом шаге
|
||||
избавляясь от выбросов.
|
||||
|
||||
\section{Задания для самостоятельного выполнения}
|
||||
\begin{enumerate}
|
||||
\item Некоторая совокупность состоит из трех групп:~$X_1$, $X_2$, и~$X_3$. Группы
|
||||
имеют следующие значения:
|
||||
\begin{verbatim}
|
||||
X1 = 35.04 35.45 35.01 34.94 34.63 35.11 34.41 35.29 35.69
|
||||
34.69 35.36 35.53 34.30 34.36 35.23
|
||||
X2 = 34.30 34.80 34.86 34.81 35.08 34.79 35.04 33.93 34.48
|
||||
34.41 33.74 34.60 34.00
|
||||
X3 = 35.17 34.21 34.78 34.65 34.16 33.62 34.53 34.12 34.82
|
||||
34.77 35.29 34.81 34.28 34.72 34.12 34.55 34.53 34.55
|
||||
\end{verbatim}
|
||||
Найдите: групповые средние (35,00, 34.53, 34.54), общее среднее (34.69), групповые
|
||||
дисперсии (0.19, 0.19, 0.16), генеральную дисперсию (0.22).
|
||||
|
||||
\item Усовершенствуйте скрипт \verb'three_s.m' так, чтобы помимо основных
|
||||
вычислений на экране отображались среднее арифметическое значение массива с
|
||||
данными, а также 95\%-й доверительный интервал по критерию Стьюдента.
|
||||
|
||||
\item Определите давление в цилиндре с газом, исходя из закона Менделеева--Клапейрона:
|
||||
$pV=mRT/\mu$, если известно, что масса газа $m=2\,$грамма, $\mu=29\,$г/моль, $R=8.31$, а объем и
|
||||
температуру газа измеряли в течение минуты, получив следующие значения:
|
||||
\begin{center}\small
|
||||
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}
|
||||
\hline
|
||||
Величина &\multicolumn{10}{|c|}{Значение}\\
|
||||
\hline
|
||||
$V$, л&2.27&2.27&2.26&2.25&2.26&2.27&2.29&2.28&2.25&2.28\\
|
||||
$T$, К&399.4&399.1&399.3&396.8&399.5&400.2&400.6&403.0&399.2&401.3\\
|
||||
\hline
|
||||
\end{tabular}
|
||||
\end{center}
|
||||
Считайте, что за это время давление газа не успело сколь-нибудь значительно измениться.
|
||||
Определите погрешности измерения величин~$V$ и~$T$. Считая, что остальные величины являются
|
||||
постоянными, определите косвенную погрешность измерения~$p$.
|
||||
|
||||
Для удобства вычислений\к создайте скрипт, позволяющий для заданного ряда данных
|
||||
получить математическое ожидание, среднеквадратичное отклонение и относительную
|
||||
ошибку\н.
|
||||
|
||||
Запишите результат в виде $p=\mean{p}\pm\sigma_p$ ($p=101\pm1\,$кПа).
|
||||
|
||||
\item Для определения емкости~$C$ неизвестного конденсатора при помощи осциллографа исследовали
|
||||
затухающий импульс, возникающий при разрядке конденсатора через резистор~$R=3\,$кОм.
|
||||
По показаниям осциллографа были записаны следующие значения тока:
|
||||
\begin{center}\small
|
||||
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|}
|
||||
\hline
|
||||
$t$, с&0.0&0.1&0.2&0.3&0.4&0.5&0.6&0.7&0.8&0.9&1.0\\
|
||||
$I$, А&1.00&0.72&0.52&0.37&0.26&0.19&0.14&0.10&0.07&0.04&0.03\\
|
||||
\hline
|
||||
\end{tabular}
|
||||
\end{center}
|
||||
Известно, что погрешность амперметра составляет
|
||||
$\sigma_I=0.01\,$А. Кроме того, известно что сопротивление резистора известно с точностью~5\%.
|
||||
Из формулы $I=I_0\exp(-t/[RC])$ определите погрешность измерения емкости конденсатора.
|
||||
|
||||
Методом наименьших квадратов определите значение емкости конденсатора, исходя из
|
||||
уравнения $t=-RC\ln I$ (97\,мкФ). Запишите ответ в виде $C=\aver{C}\pm\sigma_C$.
|
||||
|
||||
Для увеличения точности эксперимента было проведено еще одно измерение, результаты
|
||||
которого несколько отличались от предыдущих:
|
||||
\begin{center}\small
|
||||
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|}
|
||||
\hline
|
||||
$t$, с&0.0&0.1&0.2&0.3&0.4&0.5&0.6&0.7&0.8&0.9&1.0\\
|
||||
$I$, А&1.00&0.75&0.56&0.41&0.30&0.23&0.17&0.12&0.10&0.07&0.05\\
|
||||
\hline
|
||||
\end{tabular}
|
||||
\end{center}
|
||||
Проверьте нулевую гипотезу о равенстве средних в обоих опытах. Определите величину емкости
|
||||
во втором случае (112\,мкФ).
|
||||
|
||||
Столь большое различие емкостей, полученных в результате двух независимых экспериментов,
|
||||
заставило предположить, что в результате длительной эксплуатации резистор~$R$ нагрелся, что
|
||||
вызвало увеличение его сопротивления. Считая емкость конденсатора прежней, определите сопротивление
|
||||
резистора во втором случае (3.5\,кОм).
|
||||
|
||||
\item Найдите\к обоими способами\н коэффициенты $a$ и~$b$ для таблично представленной
|
||||
зависимости $y(x)$, предполагая, что она имеет линейный вид. Найдите коэффициент корреляции~$x$
|
||||
и~$y$,
|
||||
Данные представлены в таблице:
|
||||
% a=5.15, b=2.74
|
||||
% должно получиться: a=5.0644, b=3.4020
|
||||
\begin{center}
|
||||
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
|
||||
\hline
|
||||
\bf x&1&2&3&4&5&6&7&8&9&10\\
|
||||
\hline
|
||||
\bf y&7.7 & 13.7 & 22.0 & 23.1 & 23.7 & 36.7 & 35.6 & 47.8 &
|
||||
50.2 & 52.1\\
|
||||
\hline
|
||||
\end{tabular}
|
||||
\end{center}
|
||||
($a=5.1$, $b=3.4$).
|
||||
|
||||
\item Известно, что некоторая зависимость (см. таблицу ниже) имеет вид
|
||||
$y=ax\sin(x)-b\ln(x)$. Определите коэффициенты~$a$ и~$b$ и постройте
|
||||
данную кривую с более детальным отображением (на векторе \verb'[1:0.05:10]').
|
||||
Подсказка: сразу же задайте вектора~$x$ и~$y$ как столбцы; матрица~$X$ задается
|
||||
командой \verb'X=[x.*sin(x) -log(x)]'.
|
||||
\begin{center}
|
||||
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
|
||||
\hline
|
||||
\bf x&1&2&3&4&5&6&7&8&9&10\\
|
||||
\hline
|
||||
\bf y&-0.68 & 8.41 & -23.0 & -37.2 & -73.2 & -39.7 & 9.14 & 21.0 &
|
||||
7.97 & -72.5\\
|
||||
\hline
|
||||
\end{tabular}
|
||||
\end{center}
|
||||
($a=7.72$, $b=14.8$).
|
||||
|
||||
\item Составьте модель эксперимента по измерению
|
||||
амплитуды напряжения в контуре, испытывающем колебания с основной частотой
|
||||
$\Omega=1000\,$Гц и двумя гармониками $\Omega\pm\omega$, где $\omega=74\,$Гц.
|
||||
Известно, что суммарное колебание описывается приближенной
|
||||
формулой $U=a\sin(\Omega t)+b\sin(\omega t)-c\cos(\omega t)$.
|
||||
Создайте интервал времен {\tt t=[0: 0.06: 120]}. Для получения идеальных
|
||||
значений~$U$ положите $a=361$, $b=117$, $c=92$. Отношение сигнал/шум при
|
||||
получении зашумленного сигнала выберите равным~20\,дБ.
|
||||
|
||||
Восстановите значения коэффициентов~$a$, $b$ и~$c$.
|
||||
|
||||
|
||||
\end{enumerate}
|
||||
\end{document}
|
||||
38
Komp_obr_SFedU/Materials4Pract/03/2015_09_30_pf.tab
Normal file
38
Komp_obr_SFedU/Materials4Pract/03/2015_09_30_pf.tab
Normal file
@ -0,0 +1,38 @@
|
||||
Date: Sun Aug 30 02:13:37 2015
|
||||
Focus: PF
|
||||
#
|
||||
# K0 K1 K2 K3 K4 K5 K6 K7
|
||||
# A0 L k F dS Z0 d d1
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
#
|
||||
: Alpha Delta : dAlp dDel : dA dZ : A Z : Stime :
|
||||
|00:30:52.65 +53:40:12.1|+003.23 -079.0|+145.7 -080.1|-182.01 10.02|00:33:14|
|
||||
|00:32:04.58 +73:40:06.2|-000.63 -076.1|-008.8 -076.1|-180.52 30.01|00:35:46|
|
||||
|12:35:08.54 +86:20:53.5|+032.23 +077.1|-040.6 -076.7|-179.94 49.98|00:38:11|
|
||||
|12:38:12.10 +66:17:44.9|+008.12 +080.9|-052.6 -080.5|-179.73 70.01|00:40:46|
|
||||
|07:51:06.93 +45:45:49.8|-002.66 +090.0|-051.0 -081.2|-134.58 69.59|00:43:41|
|
||||
|06:05:33.51 +56:39:37.4|-007.78 +056.7|-039.4 -080.2|-134.83 49.70|00:46:03|
|
||||
|03:37:12.63 +58:36:39.2|-010.08 -010.9|-009.9 -079.4|-135.12 29.71|00:48:24|
|
||||
|01:32:34.55 +50:15:37.0|-004.94 -073.4|+159.0 -083.2|-136.94 09.62|00:51:26|
|
||||
|01:46:08.88 +42:52:51.0|-008.05 -014.0|+170.1 -085.2|-089.76 09.43|00:54:28|
|
||||
|03:28:32.99 +36:43:39.2|-006.02 +035.9|-002.6 -080.8|-089.56 29.50|00:56:53|
|
||||
|04:51:43.49 +26:07:38.8|-003.67 +071.6|-037.3 -082.2|-089.33 49.65|00:59:18|
|
||||
|06:00:16.58 +13:40:11.2|-002.14 +090.4|-049.6 -083.6|-089.56 69.49|01:01:52|
|
||||
|03:55:14.79 -14:05:23.0|-000.25 +098.3|-051.5 -085.7|-044.47 69.58|01:04:47|
|
||||
|03:16:12.86 +03:01:52.8|-001.18 +088.2|-038.6 -085.1|-044.36 49.62|01:07:14|
|
||||
|02:35:53.87 +20:01:23.9|-003.12 +072.5|-003.5 -084.8|-044.15 29.69|01:09:42|
|
||||
|01:44:26.98 +36:13:33.5|-006.39 +052.2|+184.9 -088.1|-041.82 09.57|01:12:49|
|
||||
|01:12:58.16 +33:40:55.7|-002.28 +087.6|+191.8 -085.9|+003.63 09.99|01:15:59|
|
||||
|01:15:51.16 +13:42:37.5|+000.04 +082.5|+001.6 -082.5|+001.34 29.94|01:18:35|
|
||||
|01:18:52.88 -06:26:11.3|+001.97 +082.5|-037.3 -082.8|+000.74 50.08|01:21:09|
|
||||
|23:15:10.36 +03:00:03.2|+004.54 +052.3|-040.3 -080.0|+045.74 50.29|01:29:07|
|
||||
|00:00:32.24 +20:03:29.1|+003.22 +064.7|-003.6 -079.0|+046.11 30.29|01:31:36|
|
||||
|00:56:31.33 +36:16:21.9|+002.54 +083.8|+182.1 -083.1|+047.81 10.36|01:34:34|
|
||||
|00:40:11.31 +42:45:58.2|+006.76 +045.7|+178.5 -081.3|+090.08 10.46|01:37:28|
|
||||
|23:03:17.60 +36:43:17.1|+005.69 +034.0|+002.5 -076.5|+090.40 30.39|01:39:52|
|
||||
|21:45:06.80 +26:22:15.9|+005.87 +023.6|-035.6 -077.7|+090.41 50.37|01:42:15|
|
||||
|20:41:00.09 +13:41:02.8|+006.33 +016.2|-052.6 -079.5|+090.52 70.47|01:44:54|
|
||||
|18:34:12.94 +45:46:46.2|+008.51 +017.0|-051.6 -076.5|+135.37 70.36|01:47:47|
|
||||
|20:26:11.28 +56:41:58.6|+009.79 +000.1|-038.7 -075.1|+135.17 50.26|01:50:11|
|
||||
|22:59:40.78 +58:39:16.2|+009.24 -015.8|-005.0 -074.0|+135.00 30.25|01:52:36|
|
||||
|01:08:27.11 +50:11:46.7|+008.50 -020.1|+167.4 -078.6|+133.37 10.33|01:55:28|
|
||||
93
Komp_obr_SFedU/Materials4Pract/03/getSKNcoeff.m
Normal file
93
Komp_obr_SFedU/Materials4Pract/03/getSKNcoeff.m
Normal file
@ -0,0 +1,93 @@
|
||||
function SKN = getSKNcoeff(tabname, imprefix)
|
||||
%
|
||||
% SKN = getSKNcoeff(tabname)
|
||||
%
|
||||
% Calculate SKN coefficients & plot graphs
|
||||
%
|
||||
% parameters:
|
||||
% tabname - filename with table of dA/dZ
|
||||
%
|
||||
% SKN:
|
||||
% dA = K0 + K1/tg(Z) + K2/sin(Z) - K3*sin(A)/tg(Z) + K4 *cos(delta)*cos(P)/sin(Z)
|
||||
% dZ = K5 + K6*siz(Z) + K7*cos(Z) + K3*cos(A) + K4*cos(phi)*sin(A)
|
||||
%
|
||||
% K0 = A0 - azimuth zero; K1 = L - horiz axe inclination; K2 = k - collimation error;
|
||||
% K3 = F - lattitude error of vert. axe; K4 = dS - time error
|
||||
% K5 = Z0 - zenith zero; K6 = d - tube bend; K7 = d1 - cos. tube bend
|
||||
%
|
||||
% phi = 43.6535278 - lattitude
|
||||
% t = LST - Alpha - hour angle
|
||||
% P=atan(sin(t)/(tan(phi)*cos(Del)-sin(Del)*cos(t))) - parallax angle
|
||||
%
|
||||
if(nargin == 1) imprefix = ""; endif
|
||||
[Ald Alm Als Deld Delm Dels dAl_S dDel_S dA dZ A Z STh STm STs ] = ...
|
||||
textread(tabname, "|%f:%f:%f %f:%f:%f|%f %f|%f %f|%f %f|%f:%f:%f|", ...
|
||||
60, "headerlines", 8);
|
||||
A = A*pi/180; % all angles here will be in radians
|
||||
Z = Z*pi/180;
|
||||
Al = pi*(Ald+Alm/60+Als/3600)/180; % right accession
|
||||
Delsig = Deld./abs(Deld); % declination sign
|
||||
Del = pi*Delsig.*(abs(Deld)+Delm/60+Dels/3600)/180; % declination
|
||||
phi = 43.6535278 * pi / 180; % lattitude
|
||||
t = pi*(STh+STm/60+STs/3600)/12 - Al; % hour angle
|
||||
P = atan(sin(t)./(tan(phi).*cos(Del)-sin(Del).*cos(t))); % parallax angle
|
||||
cont = 1;
|
||||
while cont
|
||||
printf("\n\n\t\t\t\tIteration %d\n\n", cont);
|
||||
onescol = ones(size(dA)); % column with ones - for less square method
|
||||
cosZ = cos(Z);
|
||||
sinZ = sin(Z);
|
||||
cosA = cos(A);
|
||||
sinA = sin(A);
|
||||
tgZ = tan(Z);
|
||||
Xmatr = [onescol sinZ cosZ cosA cos(phi).*sinA];
|
||||
K = Xmatr \ dZ;
|
||||
K5 = K(1); K6 = K(2); K7 = K(3); K3 = K(4); K4 = K(5);
|
||||
dZSKN = K5 + K6*sinZ + K7*cosZ + K3*cosA + K4*cos(phi)*sinA; % dZ by SKN
|
||||
K4fr = cos(Del).*cos(P)./sinZ; % K4 multiplier
|
||||
dASKN34 = dA + K3*sinA./tgZ - K4*K4fr; % dA components fixed by K3 & K4
|
||||
Xmatr = [onescol 1./tgZ 1./sinZ];
|
||||
K = Xmatr \ dASKN34;
|
||||
K0 = K(1); K1 = K(2); K2 = K(3);
|
||||
SKN = [K0, K1, K2, K3, K4, K5, K6, K7];
|
||||
dASKN = K0 + K1./tgZ + K2./sinZ - K3*sinA./tgZ + K4*K4fr;
|
||||
ddA = dA - dASKN;
|
||||
ddZ = dZ - dZSKN;
|
||||
sddA = std(ddA); sddZ = std(ddZ);
|
||||
mddA = median(ddA); mddZ = median(ddZ);
|
||||
printf("sigma(dda) = %f, sigma(ddZ) = %f\n", sddA, sddZ);
|
||||
%printf("mean(dda) = %f, mean(ddZ) = %f\n", mean(ddA), mean(ddZ));
|
||||
printf("median(dda) = %f, median(ddZ) = %f\n", mddA, mddZ);
|
||||
surge = find(abs(ddA - mddA) > 2*sddA);
|
||||
ssz = size(surge,1);
|
||||
if(ssz != 0)
|
||||
printf("Surges: \n")
|
||||
for i = 1:ssz
|
||||
idx = surge(i);
|
||||
printf("%f (Z = %f, A = %f)\n", ddA(idx), Z(idx)*180/pi, A(idx)*180/pi);
|
||||
endfor
|
||||
Z(surge) = []; A(surge) = []; Al(surge) = []; Del(surge) = [];
|
||||
t(surge) = []; P(surge) = []; dZ(surge) = []; dA(surge) = [];
|
||||
++cont;
|
||||
else
|
||||
cont = 0;
|
||||
endif
|
||||
endwhile
|
||||
printf("SKN coefficients: K0..K7: %f, %f, %f, %f, %f, %f, %f, %f\n", ...
|
||||
K0, K1, K2, K3, K4, K5, K6, K7);
|
||||
fg = figure;
|
||||
plot(A*180/pi, [ddA ddZ], 'o');
|
||||
legend("ddA", "ddZ");
|
||||
xlabel("A, degr"); ylabel("Remaining error: real-model");
|
||||
plotgr(sprintf("%s_%s", imprefix, "diff_vs_A"), fg);
|
||||
fg = figure;
|
||||
plot(Z*180/pi, [ddA ddZ], 'o');
|
||||
legend("ddA", "ddZ");
|
||||
xlabel("Z, degr"); ylabel("Remaining error: real-model");
|
||||
plotgr(sprintf("%s_%s", imprefix, "diff_vs_Z"), fg);
|
||||
endfunction
|
||||
|
||||
function plotgr(nm, fg)
|
||||
print(fg, '-dpdf', sprintf("%s.pdf", nm));
|
||||
print(fg, '-dpng', sprintf("%s.png", nm));
|
||||
endfunction
|
||||
31
Komp_obr_SFedU/Materials4Pract/03/three_s.m
Normal file
31
Komp_obr_SFedU/Materials4Pract/03/three_s.m
Normal file
@ -0,0 +1,31 @@
|
||||
% three_s.m
|
||||
% [ X sigma ] = three_s(x, n)
|
||||
% Производит отбор выборки x с соответствующими частотами n
|
||||
% при помощи критерия "трех сигм"
|
||||
% результат: среднее значение X и его среднеквадратичное отклонение, sigma
|
||||
|
||||
function [ X sigma ] = three_s(x, n)
|
||||
newx = []; % вспомогательный массив
|
||||
Data = [x ; n]; % совмещенный массив данных
|
||||
X = sum(x.*n)/sum(n); % среднее арифметическое
|
||||
sigma = sqrt(sum(n.*(x-X).^2)/sum(n)); % среднеквадратичное отклонение
|
||||
down = X-3*sigma; % нижняя граница доверительного интервала
|
||||
up = X+3*sigma; % верхняя граница -//-
|
||||
a = find(x < down); % a и b - массив координат, выходящих за границы
|
||||
b = find(x > up);
|
||||
while (length(a) > 0) || (length(b) > 0) % пока есть неверные значения
|
||||
Data = Data(:, find(Data(1, find(Data(1,:) >= down)) <= up)); % выбрасываем их
|
||||
x = Data(1,:);
|
||||
n = Data(2,:);
|
||||
X = sum(x.*n)/sum(n);
|
||||
for a = [1:length(n)]
|
||||
newx = [newx ones(1,n(a)).*x(a)];
|
||||
endfor
|
||||
X = median(newx);
|
||||
sigma = sqrt(sum(n.*(x-X).^2)/sum(n));
|
||||
down = X-3*sigma;
|
||||
up = X+3*sigma;
|
||||
a = find(x < down);
|
||||
b = find(x > up);
|
||||
endwhile
|
||||
endfunction
|
||||
Loading…
x
Reference in New Issue
Block a user