• tworzenie macierzy i indeksowanie elementów macierzy
• operacje macierzowe i funkcje działające na macierzach
• macierze logiczne i operacje na macierzach logicznych
• optymalizacja kodu
• wektoryzacja - zamiana pętli po elementach na mechanizmy indeksowania
• pre-alokacja - dynamiczne modyfikowanie rozmiarów macierzy jest kosztowne

Ć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

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.

• porównaj szybkość działania przy zastosowaniu prealokowanego wektora x
• przyśpiesz działanie stosując wektoryzację (pozbądź się instrukcji pętli)

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

• macierze zawierające wartości logiczne (typ logical), wartość 0 to fałsz, wartość 1 (lub niezerowa) to prawda
• indeksowanie macierzy za pomocą macierzy logicznych
• funkcje false() i true() pozwalają tworzyć macierze logiczne
• funkcja logical() zamienia macierz liczbową na logiczną

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

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

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)$

• zastąp wszystkie elementy o wartości ujemnej w macierzy X liczbą 0
• wyznacz sumę wszystkich elementów $x > 0.5$

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

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

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

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:

• sprawdź, czy zmienne zawierają wartości brakujące (oznaczone jako NaN)
• usuń wartości brakujące ze zmiennej MPG i wyznacz jej wartość średnią, odchylenie i medianę
• wyznacz wartość średnią MPG dla samochodów z roku 82 (Model_year)
• wyznacz najczęściej występujący rocznik
• ile jest różnych (unikatowych) wartości zmiennej Cylinders (zob. funkcję unique())
• wypisz w kolejności od największej do najmniejszej wartości przyśpieszenia Acceleration

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

• plot(x,y) liniowy (wykres rozrzutu)
• hist(x) histogram
• boxplot(x) skrzynkowy
• bar(x) słupkowy
• pie(x) kołowy

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)

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

plot(Acceleration, Horsepower, 'o')

## 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')

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

• narysuj histogram
• wyznacz wartość średnią $\mu$, odchylenie standardowe $\sigma$
• dodaj do wykresu z histogramem wykres prezentujący funkcję rozkładu normalnego $N(\mu, \sigma) = A \cdot \exp \left(\frac{-(x-\mu)^{2}}{2 \sigma^{2}}\right)$ (dobierz $A$ tak aby wykres jak najlepiej pokrywał się z kształtem histogramu )
• porównaj uzyskany wykres z wykresem wygenerowanym za pomocą funkcji histfit()

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

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

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