diff --git a/Komp_obr_SFedU/03_Pract.pdf b/Komp_obr_SFedU/03_Pract.pdf new file mode 100644 index 0000000..1198d65 Binary files /dev/null and b/Komp_obr_SFedU/03_Pract.pdf differ diff --git a/Komp_obr_SFedU/03_Pract.tex b/Komp_obr_SFedU/03_Pract.tex new file mode 100644 index 0000000..bf0d119 --- /dev/null +++ b/Komp_obr_SFedU/03_Pract.tex @@ -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>, следовательно, можно записать ответ: $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} diff --git a/Komp_obr_SFedU/Materials4Pract/03/2015_09_30_pf.tab b/Komp_obr_SFedU/Materials4Pract/03/2015_09_30_pf.tab new file mode 100644 index 0000000..3906eb8 --- /dev/null +++ b/Komp_obr_SFedU/Materials4Pract/03/2015_09_30_pf.tab @@ -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| diff --git a/Komp_obr_SFedU/Materials4Pract/03/getSKNcoeff.m b/Komp_obr_SFedU/Materials4Pract/03/getSKNcoeff.m new file mode 100644 index 0000000..df1a410 --- /dev/null +++ b/Komp_obr_SFedU/Materials4Pract/03/getSKNcoeff.m @@ -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 diff --git a/Komp_obr_SFedU/Materials4Pract/03/three_s.m b/Komp_obr_SFedU/Materials4Pract/03/three_s.m new file mode 100644 index 0000000..0939a39 --- /dev/null +++ b/Komp_obr_SFedU/Materials4Pract/03/three_s.m @@ -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