Spis treści.
1. Minimalne wymagania programu.
2. Pliki w pakiecie.
3. Instalacja pakietu.
4. Opis działania programu.
  4.1. Wprowadzanie danych.
    4.1.1. Podstawowe dane.
    4.1.2. Funkcje.
          Spis dostępnych funkcji elementarnych.
    4.1.3. Wpółrzędne punktów pomiarowych.
    4.1.4. Niepewności pomiarowe.
    4.1.5. Przyporządkowanie punktów poszczególnym funkcjom.
    4.1.6. Początkowe wartości poszukiwanych parametrów.
  4.2. Opis wyników i dodatkowych informacji.
    4.2.1. Wartości parametrów zadanych funkcji, wraz z niepewnościami.
    4.2.2. Macierz kowariancji parametrów.
    4.2.3. Wartości współrzędnych 'poprawionych' punktów, leżących na znalezionej krzywej.
    4.2.4. Macierz kowariancji 'poprawionych' punktów.
    4.2.5. Wykres x2=f(x1), (tylko w przypadku dwuwymiarowym).
    4.2.6. Ogólne dane o przebiegu aproksymacji.
    4.2.7. Skalowanie niepewności.
  4.3. Dokładniejszy opis metody.
Do głównej strony

1. Minimalne wymagania programu.
Program MNK.BAT wymaga MS-DOS 3.0 (lub późniejszego)
W Windows 95/98/2000/XP używać MNK.EXE

2. Pliki w pakiecie.
-readme.txt - zbiór z opisem pakietu,
-mnk.ico - plik z ikoną głównego programu dla Windows,
-MNK23-s.exe - główny program, pojedyncza precyzja obliczen (największa max. liczba punktów)
-MNK23-r.exe - główny program, rzeczywista precyzja obliczen
-MNK23-d.exe - główny program, wersja podwójna precyzja obliczen (najmniejsza max. liczba punktów)
-MNK23-demo.exe - wersja demo, rzeczywista precyzja obliczen
-mnk.bat - główny program dla systemu DOS,
-rtm.exe, dpmi16bi.ovl, egavga.bgi, litt.chr - pomocnicze programy i moduły.

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ść chibył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.
 

Na początek                                           Do głównej strony