Macierze - podsumowanie

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

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 właściwie	A ./ B	dzielenie elementów
A \ B	dzielenie macierzy (rozwiązanie )	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
randn(N, M)	macierz z wartościami losowymi z rozkładu normalnego
linspace(A, B, K)	wektor puntów na odcinku
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

Ćwiczenie Rozwiąż układ równań 1. rozwiązując równanie
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 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 wartość pola dla promienia . 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

Wektory i macierze logiczne

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

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

zastąp wszystkie elementy o wartości ujemnej w macierzy X liczbą 0

wyznacz sumę wszystkich elementów

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:

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

Wykresy statystyczne

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)

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.

narysuj histogram

wyznacz wartość średnią , odchylenie standardowe

dodaj do wykresu z histogramem wykres prezentujący funkcję rozkładu normalnego (dobierz tak aby wykres jak najlepiej pokrywał się z kształtem histogramu )

porównaj uzyskany wykres z wykresem wygenerowanym za pomocą funkcji histfit()

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 i zwróci współczynniki i określające linię prostą dopasowaną do danych i za pomocą metody najmniejszych kwadratów.

Dla danych i (gdzie ) współczynniki i wyznacz ze wzoru:

gdzie i to wartości średnie zmiennych.

Argumentami funkcji są dwa wektory i . 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