Еще о фокусе и температурах

[eddy_em]

Помимо попытки найти зависимость фокуса от температуры я еще проанализировал зависимость между температурой зеркала и подкупольного. Результаты получились из разряда "отсутствие результата — тоже результат": зависимость фокуса от температуры разных элементов крайне неявная (коэффициент корреляции мал, стандартное отклонение велико); зависимость температуры зеркала от разности температур между ним и подкупольным та же, что и была раньше (для изменения температуры зеркала на 1°C за один час в подкупольном температура должна отличаться на 10°C).


Наилучший коэффициент корреляции (-0.66) получился между фокусом и температурой зеркала. Получилось, что F = 40.22 - 0.166*T_mirror. Стандартное отклонение, правда, все равно велико: 1.49°C. Вот картинка:

Фокус и температура зеркала.

Температуры я взял из своего архива, написав утилитку для извлечения данных (на основе CGI, который используется для формирования графиков и выдачи данных в веб). Так как данные шли с разрывами (то сервер ломался, то свет отключали), я выбрал 10 наиболее длинных непрерывных кусков и проанализировал их.
Получилось, что, во-первых, почему-то корреляция между температурами зеркала и подкупольного вышла несмещенной (максимум приходится четко на 0). Понятно, что в случаях, когда включалась система охлаждения, наблюдался сдвиг максимума.

Температуры для корреляции

Корреляция температур с рисунка выше

Температуры, когда система охлаждения не включалась

Корреляция температур с рисунка выше

Ну, а для скоростей изменения температуры зеркала получались примерно одинаковые коэффициенты: около 0.01час^-1. Вот пример данных с наибольшей корреляцией:

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

А вот все (за исключением выбросов) данные:

Все данные: с 31/12/10 23:29:58 по 24/09/13 13:56:57; аппроксимация: F=0.0096*dT+0.0152

Для того, чтобы преобразовать UNIX-time и юлианские даты в обычные даты я использовал простую функцию:

function S = JDtoSTR(J)
	fmt = "%d/%m/%y %T";
	if(nargin == 2) fmt = "%d/%m/%y"; endif
	J = J(:);
	S = [];
	for i = 1:size(J,1)
		S = [S; strftime(fmt, localtime((J(i)-2440587.5)*86400))];
	endfor
endfunction

(аналогично с UNIX-time, только не надо юлианские даты в UNIX-time переводить).

А вот так меняются подписи на графиках:

Lbls=JDtoSTR(get(gca(), "xtick"));
set(gca(), "xticklabel", Lbls);

Вот так производилась "нарезка" временных интервалов на непрерывные куски:

BrokenIdxs = find([JD; 0] - [0; JD] > 1/10);
[Val, BI] = sort(diff(BrokenIdxs), 'descend');
% построение всех графиков:
for i=1:10; plot_points(i, BI, BrokenIdxs, dMir, dT, JD, MirT, DomT); endfor

А вот функция для построения графиков:

function plot_points(idx, BI, BrokenIdxs, dMir, dT, JD, MirT, DomT)
	I = BI(idx);
	first = BrokenIdxs(I);
	last = BrokenIdxs(I+1) - 1;
	iall = [first:last];
	dMir_s = dMir(iall);
	dT_s = dT(iall);
	JD_s = JD(iall);
	MirT_s = MirT(iall);
	DomT_s = DomT(iall);
	DomT_s24 = movingavg(24, DomT_s);

	Tit = sprintf("Data from %s to %s", JDtoSTR(JD_s(1)), JDtoSTR(JD_s(end)));
	plot(JD_s, [MirT_s DomT_s DomT_s24]);
	legend(["Mirror temperature"; "Indome temperature"; "Indome moving average (1 day)"]);
	xlabel("Date"); ylabel("T, \\circ{}C");
	Lbls=JDtoSTR(get(gca(), "xtick"), 1);
	set(gca(), "xticklabel", Lbls);
	title(Tit);
	print(sprintf("Temperatures_%02d.png", idx), "-dpng");
	print(sprintf("Temperatures_%02d.eps", idx), "-deps -color");
	
	corr_coeff = corr(dMir_s, dT_s)
	X=[dT_s ones(size(dT_s))];
	coeff = X\dMir_s;
	plot(dT_s, dMir_s, '.', dT, coeff(1)*dT+coeff(2));
	xlabel("T_{indome} - T_{mirror}, \\circ{}C");
	ylabel("\\Delta T_{mirror}, \\circ{}C/hour");
	sig = "";
	if(coeff(1) > 0) sig = "+"; endif
	legend(["Data"; sprintf("dT_{miror}=%f%s%f\\cdot(T_{indome}-T_{miror})\ncorrcoeff = %f", ...
		coeff(2), sig, coeff(1), corr_coeff)]);
	title(Tit);
	print(sprintf("dTM_vs_dT-%02d.png", idx), "-dpng");
	print(sprintf("dTM_vs_dT-%02d.eps", idx), "-deps -color");
	
	[Xcorr Lag] = xcorr(MirT_s, DomT_s, 'coeff');
	[Xcorr24 Lag24] = xcorr(MirT_s, DomT_s24, 'coeff');
	plot(Lag/24, [Xcorr Xcorr24])
	xlabel("Lag, days");
	ylabel("Xcorr coefficient");
	legend(["Correlation of Mirror and Indome temperatures"; "Correlation of Mirror and moving average (24h) of indome"], "location", "southeast");
	Xm = max(Xcorr);
	Shft = Lag(find(Xcorr == Xm));
	text(Shft/24, Xm, sprintf(" <-- %d hours",Shft));
	Max24 = max(Xcorr24); Shft24 = Lag(find(Xcorr24 == Max24));
	text(Shft24/24, Max24, sprintf(" <-- %d hours", Shft24));
	title(Tit);
	print(sprintf("MirDome_Xcorr_%02d.png", idx), "-dpng");
	print(sprintf("MirDome_Xcorr_%02d.eps", idx), "-deps -color");
endfunction

Ну и скользящее среднее:

function Y = movingavg(Window, Data)
% Y = movavg(Window, Data)
% moving average of *Data* by patterns of *Window* size
%
	Y = filter(ones(1,Window)/Window, 1, Data);
endfunction

Comments

[p2004r.livejournal.com]

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

[eddy-em.livejournal.com]

Да с температурами-то ничего особо интересного. Это я так, ради любопытства сделал. Интересно было проверить, можно ли однозначно вычислять новое положение фокуса, исходя из изменений температуры. Получается, что нельзя. Надо лазерный дальномер вешать.

[p2004r.livejournal.com]

скорее всего можно, но не просто из температуры :) то как крайние значения практически лежат на параллельных прямых говорит о простой линейной связи. нужен только параметр в модель --- может это в какую сторону был наклонен агрегат? есть ведь записи привязанные по времени о том как установка была ориентирована в моменты замера температуры?

[eddy-em.livejournal.com]

В принципе, да: можно еще попробовать учесть координату Z. Но так влом...