Spis treści

Macierze - podsumowanie

Operacje na macierzach

Operator Opis Operator Opis
A + B dodawanie A .+ B dodawanie
A - B odejmowanie A .- B odejmowanie
A * B iloczyn macierzy A .* B mnożenie elementów
A / B dzielenie macierzy $A B^{-1}$ (a właściwie $(B^T/A^T)^T)$ A ./ B dzielenie elementów
A \ B dzielenie macierzy $A^{-1}B $ (rozwiązanie $Ax=B$) A .\ B dzielenie elementów
A' transpozycja (sprzężona) A.' transpozycja
A ^ n potęgowanie macierzy (A * A * ... * A) A.^B potęgowanie elementami

Funkcje tworzące macierze

Funkcja Opis
zeros(N, M) macierz zerowa
ones(N, M) macierz z jedynkami
eye(N, M) macierz tożsamościowa z 1 na diagonali
rand(N, M) macierz z wartościami losowymi z rozkładu jednostajnego $[0,1]$
randn(N, M) macierz z wartościami losowymi z rozkładu normalnego $N(0,1)$
linspace(A, B, K) wektor $K$ puntów na odcinku $[A,B]$
magic(N) macierz liczb całkowitych o identycznej sumie w wierszach i kolumnach

Funkcje macierzowe

Funkcja Opis
reshape(A, N, M) zmiana wymiarów macierzy
length(A) ilość elementów wektora lub najdłuższy wymiar macierzy
size(A) ilość wierszy i kolumn macierzy
numel(A) ilość elementów
trace(A) ślad macierzy
det(A) wyznacznik macierzy
rank(A) rząd macierzy
inv(A) odwrócenie macierzy $A^{-1}$

Ćwiczenie Rozwiąż układ równań $Ax=b$ $$ \left\{\begin{array}{l} x+2 y+3 z=1 \\ 4 x+5 y+6 z=1 \\ 7 x+8 y \quad=1 \end{array}\right. $$ 1. rozwiązując równanie $x = A^{-1}b$
2. za pomocą operatora \ 

A = [1 2 3 ; 4 5 6; 7 8 0]
b = [1 ; 1; 1]
x = inv(A) * b

x1 = A \ b 
A * x1
A =

   1   2   3
   4   5   6
   7   8   0

b =

   1
   1
   1

x =

  -1.0000e+00
   1.0000e+00
  -4.1633e-17

x1 =

  -1.0000e+00
   1.0000e+00
  -3.7007e-17

ans =

   1.00000
   1.00000
   1.00000

Macierze i funkcje

Wiele funkcji Matlaba $f(x)=y$ jest przystosowana do działania na macierzach o dowolnych kształtach 

X = reshape(1:24, 4, 6)

a = sin(X)
b = 5 *  cos(X) .* sin(a) + exp(X)
X =

    1    5    9   13   17   21
    2    6   10   14   18   22
    3    7   11   15   19   23
    4    8   12   16   20   24

a =

   0.8414710  -0.9589243   0.4121185   0.4201670  -0.9613975   0.8366556
   0.9092974  -0.2794155  -0.5440211   0.9906074  -0.7509872  -0.0088513
   0.1411200   0.6569866  -0.9999902   0.6502878   0.1498772  -0.8462204
  -0.7568025   0.9893582  -0.5365729  -0.2879033   0.9129453  -0.9055784

b =

 Columns 1 through 4:

             4.73259           147.25217          8101.25916        442415.24281
             5.74721           402.10475         22028.63723       1202604.85597
            19.38931          1098.93533         59874.12309       3269015.07284
            56.84211          2980.35003        162752.63456       8886111.88011

 Columns 5 and 6:

      24154953.88173    1318815732.45003
      65659966.88446    3584912846.17585
     178482301.70134    9744803448.24378
     485165197.02438   26489122128.17479
type pole_kola2.m
x = rand(3, 4)
y = pole_kola2(x) + x  
pole_kola2.m is the user-defined function defined from: /home/marek/work/zajecia/2020.zimowy/pp2/github/pole_kola2.m

function pole = pole_kola2(r)
% Skrypt wyznacza pole koła o promieniu r

pole = pi * r.^2;

end

x =

   0.388627   0.326921   0.739871   0.325203
   0.666370   0.695561   0.803623   0.548489
   0.022186   0.087533   0.405026   0.796218

y =

   0.863104   0.662685   2.459609   0.657449
   2.061393   2.215480   2.832493   1.493604
   0.023732   0.111604   0.920391   2.787873

Ćwiczenie popraw program pole_kola.m tak aby funkcja zwracała poprawny wynik dla argumentu macierzowego. Wynikiem funkcji będzie macierz zawierająca w każdym z elementów $P_{ij}$ wartość pola dla promienia $R_{ij}$. Postaraj się nie korzystać z instrukcji pętli.

Ćwiczenie zmierz czas wykonywania poniższych instrukcji. 

n=100000
x = [];
tic
for i=1:n
    x(i) = sqrt(i);
end
toc
sum(x)

x = zeros(1, n);
tic
for i=1:n
    x(i) = sqrt(i);
end
toc
sum(x)

tic
x = sqrt([1:n]);
toc
sum(x)

n =  100000
Elapsed time is 0.445198 seconds.
ans =  21082008.97392
Elapsed time is 0.412186 seconds.
ans =  21082008.97392
Elapsed time is 0.00226307 seconds.
ans =  21082008.97392

Wektory i macierze logiczne

X = [1 , -5, 7 ; -1, 42, -12]
ind = X < 0
whos X ind
X =

    1   -5    7
   -1   42  -12

ind =

  0  1  0
  1  0  1

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        X           2x3                         48  double
        ind         2x3                          6  logical

Total is 12 elements using 54 bytes
-12# indeksowanie macierzą logiczną
a = X(ind)    
b = X( X > 0 )
X( X > 0 ) = 102
ans = -12
a =

   -1
   -5
  -12

b =

    1
   42
    7

X =

   102    -5   102
    -1   102   -12
A = true(2)
B = false(2, 3)
whos A B 
A =

  1  1
  1  1

B =

  0  0  0
  0  0  0

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        A           2x2                          4  logical
        B           2x3                          6  logical

Total is 10 elements using 10 bytes
# wartości logiczne mogą być traktowane jak liczby całkowite
A = false(2)
C = A + 1
whos A C
A =

  0  0
  0  0

C =

   1   1
   1   1

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        A           2x2                          4  logical
        C           2x2                         32  double

Total is 8 elements using 36 bytes

Ćwiczenie dlaczego poniższe wyrażenie nie działa? Spróbuj to naprawić. 

w = 1:5
w([1 0 1 0 1] == 1 ) = 0
w =

   1   2   3   4   5

w =

   0   2   0   4   0

Przydatne funkcje operujące na wartościach logicznych

Funkcja Opis
any() zwraca 1 gdy wektor posiada niezerowy element
dla macierzy zwraca wektor z wartością 1 gdy w kolumnie jest niezerowa wartość
all() zwraca 1 dla każdej kolumny której wszystkie elementy są niezerowe
find() zwraca indeksy (indeksowanie liniowe) niezerowych elementów
isequal() sprawdza czy macierze mają jednakowe elementy
isnan() sprawdza czy w macierzy są elementy NaN 
w1 = [ 0 0 1 2 ];
w2 = zeros(1,10);

a = all(w1)
b = any(w1)

c = all(w2)
d = any(w2)
a = 0
b = 1
c = 0
d = 0
E = [ 1 2 ; 0 0 ; 3 * eye(2)]'
all(E)
any(E)
E =

   1   0   3   0
   2   0   0   3

ans =

  1  0  0  0

ans =

  1  0  1  1
E
# indeksy liniowe
k = find(E)
a = E(k)
E =

   1   0   3   0
   2   0   0   3

k =

   1
   2
   5
   8

a =

   1
   2
   3
   3
E
# indeksy wierszy i kolumn
[w k] = find(E)
E =

   1   0   3   0
   2   0   0   3

w =

   1
   2
   1
   2

k =

   1
   1
   3
   4
# szukanie zerowych elementów
k2 = find(~E)
k2 =

   3
   4
   6
   7
# szukanie elementów spełniających określone warunki
k3 = find(E > 2)
k3 =

   5
   8
X = E
isequal(E, X)
X =

   1   0   3   0
   2   0   0   3

ans = 1
% to samo za pomocą funkcji all()
all(all(X == E))   
ans = 1
X(1,1) = 13

isequal(E, X)
all(all(X == E))   
~any(any(X ~= E))
X =

   13    0    3    0
    2    0    0    3

ans = 0
ans = 0
ans = 0

Operatory logiczne

Operator Opis
A < B mniejszy niż
A > B większy niż
A <= B mniejszy lub równy
A >= B większy lub równy
A == B równy sobie
A ~= B rożny od siebie
A & B operator AND elementami
A \| B operator OR elementami
~A Operator NOT 
X = magic(5)

a = X(X > 5 & X < 10)
X =

   17   24    1    8   15
   23    5    7   14   16
    4    6   13   20   22
   10   12   19   21    3
   11   18   25    2    9

a =

   6
   7
   8
   9
parzyste = X(mod(X, 2) == 0)
parzyste =

    4
   10
   24
    6
   12
   18
    8
   14
   20
    2
   16
   22
c = X(~mod(X, 2) == 0 | X < 5)
c =

   17
   23
    4
   11
    5
    1
    7
   13
   19
   25
   21
    2
   15
    3
    9

Ćwiczenie utwórz macierz X o wymiarach 10×10 posiadającą losowe wartości z rozkłady normalnego $N(0, 1)$ 

X = randn(10, 10);
X( X < 0) = 0;
sum(X(X > 0.5))
sum(X(find(X > 0.5)))
ans =  33.202
ans =  33.202

Statystyki opisowe

Funkcja Opis
min minimum (w kolumnach)
max maksimum (w kolumnach)
mean wartość średnia
std odchylenie standardowe
median mediana
mode wartość modalna (dominanta)
var wariancja
kurtosis kurtoza
skewnes skośność
sort sortowanie

min

X = randn(3,5)

a = min(X)
z = min(min(X))
X =

  -0.98080  -1.33877  -0.75984  -1.10188  -1.06108
  -1.38270   1.19402  -0.73617  -0.39035  -1.23597
   0.48693  -0.11938   0.29326   0.48180   1.14709

a =

  -1.38270  -1.33877  -0.75984  -1.10188  -1.23597

z = -1.3827
% drugi argument wyjściowy to indeksy
[a ind]= min(X)

printf('Min. wartość w kolumnie 2 to %f\n', X(ind(2), 2))
a =

  -1.38270  -1.33877  -0.75984  -1.10188  -1.23597

ind =

   2   1   1   1   2

Min. wartość w kolumnie 2 to -1.338766
% minimim w wierszach
b = min(X, [], 2)
b =

  -1.33877
  -1.38270
  -0.11938
% minimim parami elementów
c = min(X, zeros(size(X)))
c =

  -0.98080  -1.33877  -0.75984  -1.10188  -1.06108
  -1.38270   0.00000  -0.73617  -0.39035  -1.23597
   0.00000  -0.11938   0.00000   0.00000   0.00000

mean, std, var

X = randn(3,5)

% średnia wartośc w kolumnach
a = mean(X)

# średnia w wierszach
b = mean(X, 2)

# średnia ze wszystkich wartości macierzy
c = mean(X(:))
X =

   0.152208  -0.153149  -0.337226  -0.065059  -0.777124
   0.601924   0.966119  -0.381971   0.475464   1.277699
  -0.696218   0.879928   0.179555   1.982327  -1.291726

a =

   0.019305   0.564299  -0.179881   0.797577  -0.263717

b =

  -0.23607
   0.58785
   0.21077

c =  0.18752
# odchylenie standardowe
s = std(X)

# wariancja
v = var(X)
s2 = sqrt(v)   # sdt = sqrt(var)
s =

   0.65920   0.62282   0.31208   1.06102   1.35948

v =

   0.434541   0.387907   0.097396   1.125765   1.848176

s2 =

   0.65920   0.62282   0.31208   1.06102   1.35948

mode, median

X = randi(10, 5, 3)

d = median(X)
X =

   5   4   6
   8   6   6
   2   7   8
   7   6   8
   3   4   5

d =

   5   6   6
# mediana to środkowa wartość po posortowaniu
a = sort(X)
a =

   2   4   5
   3   4   6
   5   6   6
   7   6   8
   8   7   8
X = randi(3, 5, 3)
# wartość domunująca
e = mode(X)
X =

   2   1   2
   1   2   2
   1   2   2
   3   2   1
   3   1   1

e =

   1   2   2

Ćwiczenie wczytaj dane carsmall i wykonaj następujące operacje:

W Matlabie dane można wczytać polecaniem load carsmall.
Jeżeli zbiór danych nie jest dostępny (np. w Octave) to pobierz go stąd: https://www.fizyka.umk.pl/~grochu/pp2/carsmall.mat

Wykresy statystyczne

clear
load carsmall
whos
% plik z danymi (dla Octave) https://www.fizyka.umk.pl/~grochu/pp2/carsmall.mat
Variables in the current scope:

   Attr Name              Size                     Bytes  Class
   ==== ====              ====                     =====  ===== 
        Acceleration    100x1                        800  double
        Cylinders       100x1                        800  double
        Displacement    100x1                        800  double
        Horsepower      100x1                        800  double
        MPG             100x1                        800  double
        Mfg             100x13                      1300  char
        Model           100x33                      3300  char
        Model_Year      100x1                        800  double
        Origin          100x7                        700  char
        Weight          100x1                        800  double

Total is 6000 elements using 10900 bytes
hist(MPG)

png 

boxplot(MPG);
% uwaga: w octave może być wymagane odinstalowanie pakietu statystycznego
% pkg load statistics
warning: the 'boxplot' function belongs to the statistics package from Octave
Forge which you have installed but not loaded.  To load the package, run
'pkg load statistics' from the Octave prompt.

Please read <https://www.octave.org/missing.html> to learn how you can
contribute missing functionality.
error: 'boxplot' undefined near line 1 column 1
% wykres Acceleration w grupach utworzonych przez Model_Year
boxplot(Acceleration, Model_Year);   
warning: the 'boxplot' function belongs to the statistics package from Octave
Forge which you have installed but not loaded.  To load the package, run
'pkg load statistics' from the Octave prompt.

Please read <https://www.octave.org/missing.html> to learn how you can
contribute missing functionality.
error: 'boxplot' undefined near line 2 column 1

Wykres rozrzutu

plot(Acceleration, Horsepower, 'o')

png 

## Wykres rozrzutu w grupach
unique(Cylinders)
ans =

   4
   6
   8
plot(Acceleration(Cylinders == 4), Horsepower(Cylinders==4), 'o')
hold on  % zachowanie wykresu, kolejna funkcja plot() doda punty na aktyalny wykres
plot(Acceleration(Cylinders == 6), Horsepower(Cylinders==6), 'x')
plot(Acceleration(Cylinders == 8), Horsepower(Cylinders==8), 's')

png

Ćwiczenie wygeneruj wektor zawierający 1000 wartości z rozkładu normalnego o średniej 5 i odchyleniu standardowym 10.

Wykres słupkowy i kołowy

bar([sum(Cylinders == 4) sum(Cylinders == 6) sum(Cylinders == 8)])

png 

pie([sum(Cylinders == 4) sum(Cylinders == 6) sum(Cylinders == 8)])

png

Zadanie 6 Regresja liniowa

Napisz funkcję o nazwie regresja, która dla danych wektorów $x$ i $y$ zwróci współczynniki $a$ i $b$ określające linię prostą $y = a x + b$ dopasowaną do danych $x$ i $y$ za pomocą metody najmniejszych kwadratów.

Dla danych $x_i$ i $y_i$ (gdzie $i=1,2,\ldots,n$) współczynniki $a$ i $b$ wyznacz ze wzoru:

$$ a=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}} $$

$$ b=\bar{y}-a \bar{x} $$ gdzie $\bar{x}$ i $\bar{y}$ to wartości średnie zmiennych.

Argumentami funkcji są dwa wektory $x$ i $y$. Jeżeli rozmiary wektorów nie zgadzają się, funkcja zwraca stosowany komunikat błędu oraz przerywa swoje działanie. Jeżeli użytkownik poda za mało lub za dużo argumentów, wówczas funkcja również zwraca odpowiedni komunikat z błędem i kończy swoje działanie. Do wypisania komunikatów z błędem wykorzystaj funkcję error().

Przykład działania funkcji: 

x = 0:.1:10;
y = 2 * x - 10;

% dodajemy troche szumu
y = y + 3*randn(size(x));

[a b] = regresja(x, y)

plot(x, y, 'o');
hold on
plot(x, a*x+b,'-');
 
a =  2.0233
b = -10.155

png