ELIPSA50 - katalogi z danymi demo.
PROSPARA.47
HIPERB_2.50
Prosta3D
3. Instalacja pakietu.
Zawartość pakietu należy skopiować do
wybranego katalogu (nazwanego np. MNK).
Program nawet nie dotyka rejestrów
Windows.
4. Opis działania
programu.
Program MNK pozwala na wykonanie obliczeń
metodą najmniejszych kwadratów, to znaczy dopasowanie parametrów zadanej
w postaci ogólnej krzywej do zadanych punktów. Aproksymacja może być wykonywana:
- w wielu wymiarach,
- z niepewnościami pomiarowymi na każdej osi,
- z użyciem macierzy kowariancji lub współczynników korelacji,
- z możliwością użycia naraz wielu funkcji
aproksymujących,
"sklejanych" za pomocą niektórych ich
parametrów lub punktów.
Funkcje opisujące dane krzywe muszą mieć postać polegającą na sprowadzeniu wszystkich wyrażeń na jedną stronę.
Obliczane parametry powinny być oznaczone
a1,a2,a3,... ,
a współrzędne punktów - x1,x2,x3,... (zamiast np. x,y,z,...).
Poszczególne funkcje, jeżeli jest ich wiele
(patrz 4.1.1.c), mogą obejmować wspólne (niektóre) punkty pomiarowe lub
wspólne parametry (patrz 4.1.2).
W danej aproksymacji istnieją możliwości
pominięcia niektórych punktów pomiarowych, lub pominięcia którejś funkcji
(patrz 4.1.5).
W przypadku wystąpienia błędu danych,
zarówno formalnego jak i logicznego,
stosowane są komunikaty o błędach.
W przypadku użycia wejściowej macierzy kowariancji czas trwania obliczeń jest w przybliżeniu
proporcjonalny do trzeciej potęgi iloczynu liczby punktów pomiarowych i liczby wymiarów. Czas obliczeń jest w ogóle wydłużony przez nieliniowość metody, spowodowaną założeniem występowania błędów
dla wszystkich mierzonych wielkości (dla wszystkich współrzędnych kazdego punktu).
Istnieją cztery sposoby wprowadzania każdego
rodzaju danych:
a. za pomocą klawiatury;
b. korzystając z poprzednio wprowadzonych
danych w programie, jeśli ostatnio była używana aproksymacja o takich samych
podstawowych danych;
c. korzystając z danych zawartych w pliku
tekstowym w jednym z podkatalogów użytkownika; w podkatalogu tym program
umożliwia także zapisanie danych wyjściowych. Dane oddziela się spacjami.
Uwaga! Raz wybranego podkatalogu użytkownika
nie można zmienić podczas pracy z jednym zestawem podstawowych danych.
d. korzystając z danych zawartych w dowolnym
pliku tekstowym.
4.1.
Dane wejściowe wprowadzane są według następującego schematu:
4.1.1. Podstawowe
dane w nast. kolejności:
a. liczba wymiarów każdego punktu (max. 3275),
b. liczba punktów pomiarowych
(max. 8190 - wersja MNK24-s),
(max. 5460 - wersja MNK24-r),
(max. 4095 - wersja MNK24-d),
c. liczba funkcji wiążących współrzędne
poszczególnych punktów (max. 20),
d. liczba parametrów we wszystkich funkcjach
(max. 3275).
Dodatkowo musi być spełniony warunek,
że iloczyn liczby punktów i liczby wymiarów musi być mniejszy niż
16384 - wersja MNK24-s,
10922 - wersja MNK24-r,
8192 - wersja MNK24-d,
Liczby te można w czasie działania programu
modyfikować.
Nazwa odpowiedniego zbioru danych w katalogu użytkownika: basedata.dat .
4.1.2. Wprowadzenie
analitycznej postaci funkcji.
Wszystkie parametry funkcji powinny być
zapisane jako a1,a2,a3,... ,
a wszystkie współrzędne punktów - jako
x1,x2,x3,... (odpowiedniki x, y, z, ...).
POSTAĆ TA POWINNA MIEĆ FORMĘ POPRAWNEGO
WYRAŻENIA UZYSKANEGO PO SPROWADZENIU WSZYSTKICH ZMIENNYCH I STAŁYCH NA
JEDNĄ STRONĘ.
Dlatego funkcja x3=f(x1,x2) musi być wyrażona w formie uwikłanej F(x1,x2,x3)=0
bez części "=0". Oto przykład:
zastępując z = a1*x^2 + a2*y + a3 przez x3= a1*x1^2 + a2*x2 +a3
i pisząc w formie F(x1,x2,x3)=0, bez części "=0", otrzymujemy:
a1*x1^2 + a2*x2 + a3 - x3
Istnieje dzięki temu możliwość zapisania
krzywych nie będących w istocie funkcjami, nie mających łatwej do wyznaczenia
postaci, a także modyfikowania przez użytkownika takich wyrażeń funkcyjnych,
w celu uniknięcia wyrażeń nie dających się obliczyć.
Ostatni przypadek może także zajść, gdy zmienione w czasie iteracji programu
wartości współrzędnych punktów lub parametrów funkcji wykroczą poza dziedzinę
danego wyrażenia.
Funkcje powinny być ciągłe w podlegającym
obliczeniom obszarze zmienności zarówno współrzędnych punktów jak i parametrów.
Ponadto funkcje powinny być w przybliżeniu
liniowe w otoczeniu (plus-minus niepewność standardowa) każdego punktu pomiarowego i w otoczeniu
wartości każdego parametru.
Należy pamiętać, że zwłaszcza te ostatnie
wartości mogą znacznie zmienić się od początkowych do poszukiwanych w trakcie
iteracji programu .
POSZCZEGÓLNE FUNKCJE MOGĄ OBEJMOWAĆ WSPÓLNE PARAMETRY LUB WSPÓLNE PUNKTY.
Obliczenia wartości funkcji oraz potrzebnych w programie pochodnych funkcji są dokonywane z dokładnością 19 miejsc znaczących, chociaż dane są przechowywane w rzeczywistej precyzji (pojedynczej lub podwójnej, w zależności od wersji) w celu oszczędzania pamięci.
Nazwa odpowiednich zbiorów w katalogu użytkownika: function.1 , function.2 , ...
Spis działań i podstawowych
funkcji dostępnych w programie:
+,-,*,/ - działania podstawowe;
x^y - podniesienie x do potęgi y;
pi - wartość pi=3.141592653589793238 (funkcja
bezparametrowa);
sqr(x) - kwadrat argumentu;
sqrt(x) - pierwiastek kwadratowy;
exp(x) - funkcja exponencjalna;
ln(x) - logarytm naturalny;
cos(x) - cosinus (w radianach);
sin(x) - sinus (w radianach);
tg(x) - tangens (w radianach);
ctg(x) - cotangens (w radianach);
arctg(x) - arcus tangens (w radianach);
cosh(x) - cosinus hiperboliczny;
sinh(x) - sinus hiperboliczny;
tgh(x) - tangens hiperboliczny;
abs(x) - wartość bezwzględna;
sgn(x) - funkcja signum;
H(x) - skok jednostkowy (funkcja Heaviside'a);
H(x)=1 gdy x>=0, inaczej H(x)=0;
random - liczba losowa z przedziału [0,1)
lewostronnie domkniętego;
rand(n) - liczba losowa z przedziału [0,n)
lewostronnie domkniętego,
w obu przypadkach rozkład prawdopodobieństwa
jednostajny;
besj0(x) - funkcja Bessela pierwszego
rodzaju, 0-go rzędu;
besj1(x) - funkcja Bessela pierwszego
rodzaju, 1-go rzędu;
besj(n,x)- funkcje Bessela pierwszego
rodzaju, n-tego rzędu, n-naturalne;
besy0(x) - funkcja Bessela drugiego rodzaju,
0-go rzędu, x>0, (f. Webera);
besy1(x) - funkcja Bessela drugiego rodzaju,
1-go rzędu, x>0, (f. Webera);
besy(n,x)- funkcje Bessela drugiego rodzaju,
n-tego rzędu, n-naturalne, x>0, (f. Webera);
besi0(x) - zmodyfikowana funkcja
Bessela pierwszego rodzaju, 0-go rzędu, { besi0(x)=besj0(ix) };
besi1(x) - zmodyfikowana funkcja
Bessela pierwszego rodzaju, 1-go rzędu, { besi1(x)=-i*besj1(ix) };
besi(n,x)- zmodyfikowane funkcje
Bessela pierwszego rodzaju, n-tego rzędu, n-naturalne,
{ besi(n,x)=(i^n)*besj(n,ix) };
besk0(x) - zmodyfikowana funkcja
Bessela drugiego rodzaju, 0-go rzędu, x>0, (f. MacDonalda);
besk1(x) - zmodyfikowana funkcja
Bessela drugiego rodzaju, 1-go rzędu, x>0, (f. MacDonalda);
besk(n,x)- zmodyfikowane funkcje
Bessela drugiego rodzaju, n-tego rzędu, n-naturalne, x>0,
(f. MacDonalda);
gamln(x) - logarytm funkcji gamma, x>0
lub x<0 i x należy do (-2*k,-2*k+1), gdzie k naturalne;
gamm(x) - funkcja gamma, x>0 ;
gamma(x) - funkcja gamma, x rózne od -k
, gdzie k-naturalne;
gammp(a,x)-funkcja gamma niezupelna, a>0
i x>=0;
bet(x,y) - funkcja beta, x>0 i y>0 ;
beta(x,y)- funkcja beta, x,y,(x+y) różne
od -k, gdzie k-naturalne;
beti(a,b,x) -funkcja beta niezupełna,
a>0, b>0, x>=0, x<=1 ;
erf(x) - całka prawdopodobieństwa (error
function);
lgndr(n,x)-wielomiany Legendre'a, n-calkowite,
|x|<=1 ;
plgndr(n,m,x) -wielomiany Legendre'a (sferyczne
harmoniki), n,m-całkowite, |x|<=1, m>=0, m<=n ;
4.1.3. Wprowadzenie
wpółrzędnych punktów pomiarowych.
W odpowiednim pliku tekstowym wartości
współrzędnych (i tylko one) powinny być odseparowane spacją, znakiem tab lub końcem
linii i umieszczone w następującej kolejności:
dla punktu 1 : x1 x2 x3 ...
dla punktu 2 : x1 x2 x3 ... itd.
Współrzędne kolejnych punktów mogą być
wpisywane w tej samej linii.
Każda wartość powinna mieć postać liczby
rzeczywistej, której mantysa może zawierać 7 (dla trybu single), 11 (dla trybu real) lub 15 (dla trybu double) cyfr znaczących. Dla oznaczenia potęgi
liczby 10 należy używać litery E lub e.
Każda dodatnia wartość powinna być większa niż
ok. 2.9E-39 i mniejsza niż 1.7E38 (odpowiednio 5.0E-324 i 1.7E308 dla podwójnej precyzji).
Jednak ze względu na obliczenia wartości
wyrażeń dla wprowadzonych funkcji oraz ograniczenie z punktu 4.1.4, zakres
ten prawdopodobnie będzie mniejszy. Dlatego wtedy pożądane jest takie przeformułowanie
problemu, żeby wartości te mieściły się bliżej jedności.
Nazwa odpowiedniego zbioru w katalogu użytkownika: points.dat .
4.1.4. Wprowadzenie
niepewności pomiarowych.
Założeniem metody, na której opiera się
program, jest podleganie wielkości mierzonej (z błędem) rozkładowi normalnemu
Gaussa.
Postać danych powinna być taka, jak w
punkcie 4.1.3 .
Istnieją dwie możliwości wprowadzenia
niepewności :
1. Odchyłki standardowe
dla każdego wymiaru każdego punktu, przy czym za pomocą klawiatury można
ten proces skrócić, jeżeli odchyłki są takie same dla jednego wymiaru wszystkich
punktów.
Ze względu na obliczanie przez program
kwadratów odchyłek, ich dopuszczalny zakres wynosi około: od 1.E-19 do 1.E19 (odpowiednio 1.E-162 i 1.E154 dla podwójnej precyzji).
Kolejność wprowadzania niepewności pomiarowych
- taka, jak w punkcie 4.1.3 .
Nazwa odpowiedniego zbioru w katalogu użytkownika:
errors.dat .
Współczynniki korelacji mogą być uwzględnione dla różnych wymiarów każdego punktu.
Mogą one być użyte, gdy są skorelowane błędy różnych wymiarów dla każdego punktu osobno.
Kolejność wprowadzania dla każdego punktu w czterowymiarowym przypadku jest następująca:
r12, r13, r14, r23, r24, r34.
Nazwa odpowiedniego zbioru w katalogu użytkownika: correls.dat .
2. Macierz kowariancji punktów pomiarowych
(tylko z pliku tekstowego).
Może być użyta zwłaszcza wtedy, gdy skorelowane są błędy różnych punktów.
Jej rozmiar wynosi n*n (gdzie n jest iloczynem
liczby punktów i liczby wymiarów), a główna przekątna składa się z kwadratów
odchyłek standardowych.
Kolejność musi odpowiadać kolejności z
punktu 4.1.3 .
Nazwa odpowiedniego zbioru w katalogu użytkownika:
covars.dat .
4.1.5. Wprowadzenie
przyporządkowania punktów pomiarowych poszczególnym funkcjom (krzywym).
Nie ma takiej konieczności, gdy w programie
użyta jest jedna funkcja.
Postać danych w odpowiednim pliku tekstowym
powinna być następująca:
w każdym wierszu umieszczone numery punktów
przyporządkowanych odpowiedniej funkcji, przy czym kolejność wierszy taka,
jak podczas wprowadzania funkcji, n.p. (jeśli trzy funkcje są użyte):
1-10 20 21 22
0
5 7-12
Istnieją możliwości pominięcia niektórych
punktów pomiarowych, uwzględnienia danego punktu kilka razy w różnych funkcjach,
lub pominięcia którejś funkcji przez umieszczenie zera w odpowiadającym
jej wierszu .
Nazwa odpowiedniego zbioru w katalogu
użytkownika: whichpnt.dat .
4.1.6. Wprowadzenie
początkowych wartości poszukiwanych parametrów.
Użyta w programie metoda dla rozpoczęcia
swojej iteracji wymaga podania początkowych wartości poszukiwanych parametrów.
Jeżeli przyjęta funkcja znacznie różni się od funkcji liniowej, to iteracja
będzie zbieżna jedynie wtedy, jeżeli wartości początkowe są dobrymi przybliżeniami
wartości rzeczywistych. Nie można podać tu żadnej ogólnej reguły wyboru
tych wartości, bowiem zależy ona od natury problemu.
Wartości te powinny mieć postać liczby
rzeczywistej z ograniczeniami jak w punkcie 4.1.3.
Kolejność wprowadzenia musi odpowiadać
kolejności oznaczenia parametrów przyjętej podczas wprowadzenia funkcji.
Nazwa odpowiedniego zbioru w katalogu
użytkownika: params.dat .
4.2. Obliczone dane
wyprowadzane są według następującego schematu:
4.2.1. Wartości parametrów
zadanych funkcji, wraz z niepewnościami
(na ekran i do pliku aprxpars.dat ) .
4.2.2. Macierz kowariancji
parametrów
(gdy jest żądana, lecz tylko do pliku
parcovar.dat ).
4.2.3. Wartości
współrzędnych 'poprawionych' punktów, leżących na znalezionej
krzywej
(na ekran i do pliku aprxpnts.dat ).
4.2.4. Macierz kowariancji
'poprawionych' punktów
(gdy jest żądana, lecz tylko do pliku
pntcovar.dat ).
4.2.5. Wykres x2=f(x1)
, (tylko w przypadku dwuwymiarowym).
Istnieje możliwość rysowania prostokątów
lub elips niepewności oraz wykresu znalezionej krzywej, traktując ją jako funkcję
jednej ze zmiennych x1,x2 lub jako dowolną krzywą uwikłaną. W tym ostatnim
przypadku czasami wykres może być niezadowalający i wtedy poprawa może
nastąpić po poszerzeniu zakresu danej krzywej. Możliwość dobrania zakresu
dla wykresu każdej linii ma zapobiegać błędowi wynikłemu z przekroczenia
dziedziny funkcji. Zapisanie wykresu w pliku bmp za pomocą skrótu Ctrl S.
4.2.6. Ogólne dane
o przebiegu aproksymacji.
- Wartość chi2 (o rozkładzie
chi-kwadrat);
powinna być równa (w granicach podanej
niżej odchyłki) liczbie stopni swobody (patrz niżej).
Za duża wartość tej funkcji świadczy o
niewłaściwym dobraniu analitycznej postaci krzywej do punktów pomiarowych,
lub o przyjęciu za małych niepewności pomiarowych. Za mała wartość świadczy
o założeniu zbyt dużej liczby parametrów, lub o przyjęciu za dużych niepewności
pomiarowych.
- Liczba stopni swobody; równa liczbie
punktów pomiarowych pomniejszonej o liczbę poszukiwanych parametrów. Liczba
ta jest też wartością oczekiwaną rozkładu wartości chi2.
- Odchyłka standardowa rozkładu wartości
chi2.
- Prawdopodobieństwo,
że, teoretycznie, wartość chi2 byłaby większa (mniejsza) niż jej wartość
otrzymana w danej aproksymacji metodą najmniejszych kwadratów.
- Współczynnik korelacji "ta-sama-strona-krzywej", razem z jego odchyłką standardową. Gdy punkty pomiarowe położone są przypadkowo po obu stronach krzywej (powierzchni, hiper-powierzchni), to współczynnik ten przyjmuje wartość bliską zeru. Gdy punkty układają się w grupy, z których każda położona jest po jednej stronie krzywej, to współczynnik znacznie odbiega od zera, t.zn. jest większy niż jego podwójna odchyłka standardowa, lub nawet zbliża się do jedności, wskazując, że krzywa została niewłaściwie wybrana (może potrzebnych jest więcej parametrów) albo punkty są skorelowane (może punkty z jednej grupy pochodzą z jednego miernika z biasem). Gdy współczynnik jest ujemny i znacznie odbiega od zera lub nawet zbliża się do minus jedności, jest to statystycznie mało prawdopodobne (może punkty zostały stworzone "ludzką ręką").
- Numery punktów, dla których różnica
jakiejkolwiek współrzędnej punktu pomiarowego i odpowiedniej współrzędnej
punktu poprawionego jest większa niż odpowiednia potrójna niepewność pomiarowa.
4.2.7. Skalowanie niepewności.
Jeżeli jesteś pewny dopasowywanej krzywej i stosunków niepewności dla wszystkich punktów i wymiarów, zamiast samych niepewności, możesz pomnożyć (menu Opt) wszystkie wejściowe i wyjściowe niepewności (punktów i parametrów) przez pierwiastek ze stosunku wyjściowego parametru chi2 do liczby stopni swobody. W ten sposób możesz uzyskać najlepsze estymaty niepewności w tym przypadku niepełnej informacji o wejściowych niepewnościach.
4.3. Dokładniejszy
opis metody.
Metoda polega na minimalizacji wielkości
M, równej
M = dT*Gy*d
gdzie znak T oznacza
transponowanie, d jest wektorem różnic współrzędnych
[ (X1-Xo1) , (Y1-Yo1) , (X2-Xo2)
, (Y2-Yo2) , ... ] ,
Xoi,Yoi oznaczają
współrzędne i-tego punktu pomiarowego,
Xi,Yi - współrzędne
i-tego "poprawionego" punktu znajdującego się na krzywej,
a Gy jest odwrotnością macierzy
kowariancji (patrz p. 4.1.4).
W przypadku dwuwymiarowym i z zerowymi
kowariancjami wielkość M jest równa
_____
\
\
(Xi-Xoi)2
(Yi-Yoi)2
/
---------- + ----------
/_____ (Sxi)2
(Syi)2
i
gdzie sumowanie odbywa się po numerach
punktów i,
a Sxi,Syi
- odpowiednie odchyłki standardowe i-tego punktu.
Wartość chi2 (z punktu 4.2.6)
jest równa (w przypadku dwuwymiarowym i z zerowymi kowariancjami)
_____
\
\
[ -F'(Xi)*(Xi-Xoi)
+ (Yi-Yoi) ]2
/
--------------------------------------
/_____
[F'(Xi)]2 * (Sxi)2
+ (Syi)2
i
gdzie F'(Xi) oznacza
cząstkową pochodną po x szukanej funkcji,
liczoną w i-tym punkcie; reszta oznaczeń
- jak wyżej.
Iteracja jest przerywana, gdy różnica między dwiema
wartościami M z kolejnych iteracji jest mniejsza niż 0.1*sqrt(M) oraz gdy każda różnica między dwiema
wartościami każdego parametru wziętymi z kolejnych iteracji jest mniejsza niż odpowiednia odchyłka standardowa parametru pomniejszona 10 razy.