### Vyhledávání

library(CRANsearcher)
#CRANsearcher()

library(pkgsearch)
#pkg_search_addin(query = "", viewer = c("dialog", "browser"))


### Základní matematické funkce

x <- -7
x = -7

# absolutni hodnota (absolute value)
abs(x)

# nasobeni a deleni (multiplication and division)
3*3
x*3
9/3
x/3
(x-3)/(x/3)

1/0
log(0)
Inf*Inf
Inf/Inf
0/0

# celociselne deleni (integer division)
13%/%2
library(numbers)
div(13,2)  
# zbytek po deleni (reminder after division)
13%%2
mod(13,2) 

# Prevod desetinných císel na zlomky (Convert decimal numbers to fractions) ###

library(MASS)
fractions(-0.1666667)
fractions(0.14)
fractions(0.4)

library(schoolmath)
decimal2fraction(0.1,6)

library(fractional)
fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE)
numerators(-0.1666667) # vypise citatele (numerator)
denominators(-0.1666667) # vypise jmenovatele (denominator)
numerical(fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE))
unfractional(fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE))

# zjednodusení zlomku (fractions simplification) ###
library(schoolmath)
cancel.fraction(42, 56)
cancel.fraction(-167, 100)


# mocniny
x = 9
y=4

sqrt(x)
x^(1/2)
x^(-1/2)

x**2
x^2
x^(-2)

library(numbers) # Determine whether p is the power of an integer.
isIntpower(144)
isIntpower(8)


# goniometricke funkce
cos(x)
sin(x)
tan(x)
acos(x)
asin(x)
atan(x)
atan2(y, x)
cospi(x) # cos(pi*x)
sinpi(x) # sin(pi*x)
tanpi(x) # tan(pi*x)

# logaritmy
x = 10
log(x)  # base e
log(x, base = exp(1))  # base e
logb(x, base = exp(1)) # base e
log10(x) # base 10
log(x,base = 10)
log2(x)  # base 2
log(x,base = 2)
log1p(x) # computes log(1+x) 


10^x # exponent
10**x # exponent
exp(x)   # e^x


library(schoolmath)
x <- c(1.002, 2, 3.14, 4, 5, 6.31, 7) 
schoolmath::is.decimal(x)
schoolmath::is.whole(x)

is.integer(x)
library(ttutils)
isInteger(x)
isInteger(c("test", 1, 2, 2.1)) # vektor (vector)

library(litteR)
is_natural_number(12)
is_natural_number(1.2)

library(numbers)
numbers::isNatural(x)


x <- c(1, 2, 3, 4, 5, 6, 7) 
schoolmath::is.even(x)
schoolmath::is.odd(x)
x %% 2 == 0 
x %% 2 != 0
is.even2 <- function(x) x %% 2 == 0; is.even2(x) 
is.odd2 <- function(x) x %% 2 != 0; is.odd2(x) 

x <- c(-1, -2, 3.02, 4, -5.2, 6, -7, 0) 
schoolmath::is.negative(x)
schoolmath::is.positive(x)
schoolmath::is.real.positive(x)
x>0
x<0
x==0


# Prvocisla
x <- c(1, 2, 3, 4, 5, 6, 7) 
schoolmath::is.prim(x)
library(DescTools)
DescTools::IsPrime(x)
library(numbers)
numbers::isPrime(x)

### Nejvetsí spolecný delitel (greatest common divisor) ###
x = 12
y = 8

library(schoolmath)
schoolmath::gcd(x, y)

library(FRACTION)
FRACTION::gcd(x, y)

library(DescTools)
DescTools::GCD(x, y, na.rm = FALSE)
DescTools::GCD(c(10,12,8), na.rm = FALSE)

library(numbers)
numbers::GCD(12, 10) 
numbers::GCD(46368, 75025)
numbers::mGCD(c(2, 3, 5, 7) * 11)
numbers::mGCD(c(2*3, 3*5, 5*7))

numbers::coprime(12, 10) # maji/nemaji prvociselneho delitele
numbers::coprime(46368, 75025) 


### Nejmensí spolecný násobek (smallest common multiple) ###
x = 12
y = 8

library(schoolmath)
schoolmath::scm(x, y)

library(DescTools)
DescTools::LCM(x, y)
DescTools::LCM(c(12,10,8), na.rm = FALSE)

library(numbers)
numbers::LCM(12, 10) 
numbers::LCM(46368, 75025) 
numbers::mLCM(c(2, 3, 5, 7) * 11) 
numbers::mLCM(c(2*3, 3*5, 5*7))


### Vektory


1:6
c(2:12)
rev(c(2:12))
seq(from=1, to=20, by=2)
seq(5, -5, -1)
seq(from=5, by=-1, along.with = 1:20)
rep(1, len = 9)
rep(1:4, len = 9)
rep(1:4, 2)
rep(1:4, each = 2)       # not the same.
rep(1:4, c(2,2,2,2))     # same as second.
rep(1:4, c(2,1,2,1))
rep(1:4, each = 2, len = 4)    # first 4 only.
rep(1:4, each = 2, len = 10)   # 8 integers plus two recycled 1's.
rep(1:4, each = 2, times = 3)  


sample(12:99, 9, replace = FALSE) # bez opakovani
sort(sample(1:20, 9, replace = FALSE))
sample(12:99, 9, replace = TRUE) # s opakovanim
sort(sample(1:20, 9, replace = TRUE))

set.seed(123)
x = sample(12:99,19, replace = FALSE);x
is.vector(x)
x
x[1]
x[1:5]
x[c(1:5)]
x[length(x)]
x[(length(x)-1)]
x[-1]
x[-(1:5)]
x[-c(1:5)]
x[-length(x)]
x[length(x)-3:length(x)]#poslednich 6 hodnot
x[c(1,4,7)]#vybrane hodnoty
which.max(x)
which.min(x)
library(Rmpfr)
Rmpfr::which.min(x)
which(x<=50)
x[which(x<=50)]
x[x<=50]
which(x>=50)
x[which(x>=50)]
which(x == max(x))#maximalni hodnota
which(x == 10)
which(x != 10)
x[x > 40]#hodnoty vyssi nez 400.
x[x < 40]#hodnoty mensi nez 100
which(x > 10)#hodnoty vyssi nez 10.
which(x < 10)#hodnoty mensi nez 10.
seq_along(x)

rank(x)  # vhodne pro vektory, nevhodne pro matice 
is.unsorted(x)
sort(x)
sort(x, decreasing = FALSE)
sort(x, decreasing = TRUE)
sort(x[which(x>=50)])
order(x)  # poradi (indexy) pro sort
x[order(x)]

a = c(1, 3, 5, 7) 
b = c(1, 2, 4, 8) 
5 * a
a + b 
a - b 
a * b 
a / b

u = c(10, 20, 30) 
v = c(1, 2, 3, 4, 5, 6, 7, 8, 9) 
u + v 

max(x)
min(x)
range(x)
library(DescTools)
Range(x)
sum(x)
cumsum(x)
prod(x)
cumprod(x)
diff(x)
sign(x)


### Matice

x1 = sample(12:99,19, replace = FALSE)
x1
x2 = sample(12:99,19, replace = TRUE)
x2
x12 = rbind(x1,x2)
x12 = cbind(x1,x2)
is.matrix(x12)
nrow(x12)
ncol(x12)

x12[order(x1),]
x12[order(x2),]
x12[,1]
x12[,"x2"]

colnames(x12)
colnames(x12)=c("X","Y")
rownames(x12)
rownames(x12) = LETTERS[1:length(x1)] 
t(x12)


### Dataframes


x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x12 = rbind(x1,x2,z)
x12 = cbind(x1,x2,z)
nrow(x12)
ncol(x12)
is.data.frame(x12)
is.matrix(x12)

x12 = as.data.frame(x12)
is.data.frame(x12)
colnames(x12)
rownames(x12)
as.matrix(x12)


mean(x1)
mean(as.numeric(x12[,1]))


as.character(x1)

m = cbind(x1,x2,x1,x2)
apply(m,1,mean)


### Seznamy (lists)

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x123 = list(x1,x2,x3,z)
is.list(x123)
names(x123) = c("x1","x2","x3","z") 
x123$x1
x123[[1]]


x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)
LL = list(x1,x2,x3)
names(LL) = c("x1","x2","x3")

lapply(LL, mean)
sapply(LL, mean)


unlist(LL)
do.call(rbind,LL)
do.call(cbind,LL)

stack(LL)


### Apply funkce
# https://www.r-bloggers.com/r-tutorial-on-the-apply-family-of-functions/
# http://www.datasciencemadesimple.com/apply-function-r/

### lapply a sapply (vektory, seznamy) ###

### Vypocitejte molekulove hmotnosti alkanu C1 - C12.
# lapply
Mhc = lapply(c(1:12), function(n){n*12.011 + (2*n + 2)*1.008})
Mhc = unlist(Mhc)
names(Mhc) = c(1:12)
Mhc
# sapply
Mhc = sapply(c(1:12), function(n){n*12.011 + (2*n + 2)*1.008})
names(Mhc) = c(1:12)
Mhc

library(live)
data(wine)
wine
head(wine)

# mean pomoci sapply a lapply
wn = sapply(c(1:length(wine[1,])), function(x){mean(wine[,x])})
names(wn) = colnames(wine)
wn
#
wn = lapply(c(1:length(wine[1,])), function(x){mean(wine[,x])})
names(wn) = colnames(wine)
wn

### tapply (matice, dataframes) ###

tapply(wine$alcohol,wine$quality,mean)
tapply(wine$alcohol,wine$quality,median)
tapply(wine$alcohol,wine$quality,length)
tapply(wine$alcohol,wine$quality,boxplot.stats)

### apply (matice, dataframes) ###

apply(wine,1,length)
apply(wine,2,length)
apply(wine,2,mean)
apply(wine,2,median)


### Matematické a fyzikální konstanty


### Ludolfovo cislo ####
pi
print(pi, digits=20)

### Eulerovo cislo ####
exp(1)
print(exp(1), digits=20) 

library(marelac)
data(Constants)
Constants

Constants$gasCt1
Constants$gasCt2
Constants$gasCt3
R = as.numeric(Constants$gasCt2[1]); R # hodnota
Constants$gasCt2[2] # jednotka


library(constants)
codata
lookup("planck", ignore.case=TRUE)
lookup("avogadro", ignore.case=TRUE)
lookup("faraday", ignore.case=TRUE)
lookup("molar", ignore.case=TRUE)
lookup("molar mass", ignore.case=TRUE)
lookup("light", ignore.case=TRUE)
lookup("boltzmann", ignore.case=TRUE)
lookup("mass", ignore.case=TRUE)
lookup("mass constant", ignore.case=TRUE)
lookup("proton mass", ignore.case=TRUE)
lookup("electron mass", ignore.case=TRUE)

# Univerzalni plynova konstanta
lookup("gas", ignore.case=TRUE)
codata[134,]
R = as.numeric(codata[134,5]) # hodnota
Ru = as.numeric(codata[134,6]) # nejistota
codata[134,7] # jednotka



### Prevody jednotek

library(scales) # %
percent(0.05)

library(pracma) # uhly vs radiany
deg2rad(120) 
rad2deg(3.14)


library(measurements)
conv_unit_options # pouzivane jednotky

# length
conv_unit(1, from="mm", to="m")
conv_unit(1, from="nm", to="angstrom")

# pressure
conv_unit(101, from="Pa", to="mmHg")
conv_unit(1, from="atm", to="Pa")

# temperature
conv_unit(100, from="C", to="K")

# energy
conv_unit(500, from="cal", to="J")
conv_unit(1, from="erg", to="J")

#volume
conv_unit(100, from="l", to="m3")

# mass
conv_unit(100, from="g", to="mg")
conv_unit(100, from="ug", to="g")

# speed
conv_unit(1, from="km_per_hr", to="m_per_sec")
conv_unit(1, from="light", to="km_per_hr")
conv_unit(1, from="light", to="m_per_sec")


### multiunits
conv_multiunit(x = 1, from="ug / l", to="g / m3")
conv_multiunit(x = 1, from="mg / l", to="kg / m3")
conv_multiunit(x = 1, from="cal / kg", to="J / g")


#### nasobky jednotek SI
library(sitools)
f2si(10000)
f2si(0.023, unit="l")
f2si(3.5e-5, unit="mol")
numbers <- c(1e5, 3.5e-12, 0.004) 
f2si(numbers, unit="g")

sapply(numbers,function(x){f2si(x, unit="g")})


# Prevod casovych jednotek

library(datetime)
xx = as.minute(60) # 60 min
as.second(xx)
as.hour(xx)
as.day(xx)

dur = duration(week=0,day=0.04166667,hours=1,minutes=0,seconds=3600)
as.numeric(dur, "minutes")


### Plynný systém mení svuj objem o 1200 ml za konstantního vnejsího tlaku 30 atm.Jakou práci vykoná plyn pri expanzi? Vyjádrete v ruzných energetických jednotkách.
# W = p . dV
dV = 1200 # [ml]
dV = conv_unit(dV, from="ml", to="m3"); dV
p = 30 # [atm] to [Pa]
p = conv_unit(p, from="atm", to="Pa"); p
W = p*dV; W # [J]
W1 = conv_unit(W, from="J", to="cal"); W1 # [cal]
W2 = conv_unit(W, from="J", to="Wsec"); W2 # [W.s]
W3 = conv_unit(W, from="J", to="erg"); W3 # [erg]


### Kolik tepla se uvolní pri pruchodu 26.43 coulombu elektrického proudu vodicem pri potenciálovém spádu 2.432 V.
# 64.33 J, 64.33e7 erg, 15.37 cal
# Q = U . I . t = U . q 
q = 26.43 # [C]
U = 2.432 # [V]
Q = U*q; Q # [J]
Q1 = conv_unit(Q, from="J", to="cal"); Q1 # [cal]
Q3 = conv_unit(Q, from="J", to="erg"); Q3 # [erg]


### Pokus trval 7658 minut. Kolik je to dní, hodin a minut?
tt = 7658 # [min]

library(lubridate)
seconds_to_period(tt*60)
library(datetime)
seconds_to_period(as.second(as.minute(tt)))


##### Prvky

library(PeriodicTable)
data(periodicTable)
periodicTable

library(CHNOSZ)
# atomova hmotnost
AW = mass("Ca"); AW
# molarni hmotnost
MW = mass("CaCO3"); MW
# pocet atomu prvku ve vzorci
ce = count.elements("CaCO3"); ce
as.chemical.formula(ce, drop.zero = TRUE)

library(IsoSpecR)
data(isotopicData)
isotopicData

library(loon.data)
data(elements)
elements


#### Zaokrouhlování na daný pocet desetinných míst

round(1234.567,2)
round(123.456,digits=2)
ceiling(1234.567)
floor(1234.567)

library(guf)
round_something(123.456, decimals = 2)
round_something("123.456", decimals = 2)

library(PKNCA)
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "e")
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "x10^")
roundString(3141.59265, digits = 3, sci_range = Inf, sci_sep = "e")


### Zaokrouhlování na daný pocet platných císlic

signif(1234.567,4)
signif(123.456,digits=4)

library(PKNCA)
signifString(3141.59265, digits = 1, sci_range = 0, sci_sep = "e")
signifString(3141.59265, digits = 1, sci_range = 0,sci_sep = "x10^")
signifString(3141.59265, digits = 1, sci_range = Inf, sci_sep = "e")


### NA and data imputation 1

### vektory

library(datasets)
head(airquality)
summary(airquality)
airquality[,"Ozone"]
airquality[,"Solar.R"]

x = airquality[,"Ozone"]
length(x)
any(is.na(x))
is.na(x[5])
is.na(x[6])

xx[6]<- NA
xx[6]

library(jwutil) # rel. podil hodnot ve vektoru
propIsNa(x)

xx = na.omit(x)
length(xx)

xx[6]<- NA
xx[6]


### Matice a dataframes

x = airquality
any(is.na(x))
is.na(x)

library(splus2R)
splus2R::which.na(x)

library(GLDEX)
GLDEX::which.na(x)

apply(x,2,splus2R::which.na)


library(ryouready)
count_na(x, along = 2)
# along: 1 = rows, 2=columns

xl = apply(x,2,splus2R::which.na)
sapply(xl,length)

library(jwutil) # rel. podil hodnot ve sloupcich
propNaPerField(x)

library(na.tools)
all_na(x)
any_na(x)
is_na(x)
which_na(x[,1])
which_na(x[1,])


library(agricolae)
delete.na(x, alternative="less")
delete.na(x, alternative="greater")
# "less"    = vynechava sloupce
# "greater" = vynechava radky

xx = delete.na(x, alternative="greater")
xx
xx[5,5]<- NA


library(jwutil) # NA na 0
zero_na(x)
zero_na(x, na_ish = TRUE)

library(DescTools)
ZeroIfNA(x)
NAIfZero(x)
#
BlankIfNA(x, blank="")
BlankIfNA(x, blank="UDL")
NAIfBlank(x)

# nahrazeni medianem sloupce
library(DescTools) 
x1 = Impute(x[,1], FUN = function(y) mean(y, na.rm = TRUE))
x2 = Impute(x[,1], FUN = function(y) median(y, na.rm = TRUE))
cbind(x[,1],x1,x2)

x[,1] = x1 # nahrazeni

library(na.tools)
na.constant(x[,1], 5)
na.max(x[,1])
na.min(x[,1])
na.mean(x[,1])
na.median(x[,1])
na.quantile(x[,1], prob=0.4)
na.mode(x[,1])
na.most_freq(x[,1])
na.inf(x[,1])
na.neginf(x[,1])
na.true(x[,1])
na.false(x[,1])
na.zero(x[,1])

set.seed(123)
xb = na.bootstrap(x[,5])
xr = na.resample(x[,5])
rbind(x[,5],xb,xr)


### opakujici se hodnoty (ties)

duplicated(c("a","b","c", "a"), incomparables = FALSE)

library(guf)
count(c("a","b","c", "a"), group.by = NULL, split.by = NULL, row.total = FALSE)


### Nejistoty a chyby mereni

### Zapis chyby (nejistoty)
library(errors)
xe <- set_errors(3.602, 0.008)
options(errors.notation = "plus-minus") 
print(xe, digits = 2)
options(errors.notation = "parenthesis")
print(xe, digits = 2)

# elementární náboj
e <- set_errors(1.6021766208e-19, 0.0000000098e-19) 
options(errors.notation = "plus-minus") 
print(e, digits = 2)
options(errors.notation = "parenthesis")
print(e, digits = 3)

# Ze zásilky kaprolaktamu bylo odebráno 10 vzorku a byl stanoven jejich bod tání. Vypocítejte hodnotu bodu tání v zásilce a její smerodatnou odchylku.
xc = c(68.5, 68.7, 68.3, 68.8, 68.5, 68.2, 68.6, 68.4, 68.2, 68.7)
mean(xc)
sd(xc) # sd
sd(xc)/sqrt(length(xc)) # str chyba prumeru
e <- set_errors(mean(xc), sd(xc)/sqrt(length(xc))) 
options(errors.notation = "plus-minus"); print(e, digits = 2)
options(errors.notation = "parenthesis"); print(e, digits = 2)

### Teplota vody byla merena dvema teplomery s ruznou presností
t1 = c(20.8, 21.1, 20.9)
t2 = c(20.95, 20.92, 20.94)
e1 <- set_errors(mean(t1), sd(t1)) 
options(errors.notation = "plus-minus"); print(e1, digits = 1)
options(errors.notation = "parenthesis"); print(e1, digits = 1)
round(100*sd(t1)/mean(t1),3) # relativni chyba [%]
e2 <- set_errors(mean(t2), sd(t2)) 
options(errors.notation = "plus-minus"); print(e2, digits = 1)
options(errors.notation = "parenthesis"); print(e2, digits = 1)
round(100*sd(t2)/mean(t2),3) # relativni chyba [%]


#### Sireni chyb ###################################

##### mean + S.E. 
## Error propagation using Taylor expansion of first- (Mean.1, sd.1) and second-order (Mean.2, sd.2)
## Monte Carlo (simulace)

## Priklad (bez uvedení poctu stupnu volnosti, resp. poctu opakovani)
library(propagate)
EXPR1 <- expression(x/y) 
x <- c(5, 0.01) 
y <- c(1, 0.01) 
DF1 <- cbind(x, y)
RES1 <- propagate(expr = EXPR1, data = DF1, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000) 
RES1
RES1$data
RES1$prop
RES1$sim
summary(RES1)
plot(RES1)


### Vypocítejte hustotu a její chybu pro látku, u níz byla opakovaným merením stanovena hmotnost 6.824 (0.008) g a objem 3.03 (0.01) ml.
EXPR2 <- expression(m/V) 
m = c(6.824, 0.008) #[g] 
V = c(3.03, 0.01) #[ml]
DF2 <- cbind(m, V) 
RES2 <- propagate(expr = EXPR2, data = DF2, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES2
RES2$data
RES2$prop
RES2$sim
library(errors)
e <- set_errors(RES2$sim[1], RES2$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)


## Priklad (s uvedením poctu stupnu volnosti) 
EXPR2 <- expression(x/y) 
x <- c(5, 0.01, 12) 
y <- c(1, 0.01, 5) 
DF2 <- cbind(x, y) 
RES2 <- propagate(expr = EXPR2, data = DF2, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES2
RES2$data
RES2$prop
RES2$sim
summary(RES2)
plot(RES2)


### Na píst o prumeru (200 +- 0.05) mm pusobí pára tlakem (8.2 +- 0.1) atm. Jakou silou pusobí pára na píst? Urcete také relativní chybu výsledku.
library(measurements)
d = c(200, 0.05) # [mm] na [m]
d = as.vector(conv_unit(c(200, 0.05), from="mm", to="m"))
p = c(8.2, 0.1)  # [atm] na [Pa]
p = as.vector(conv_unit(c(8.2, 0.1), from="atm", to="Pa"))
pi
EXPR7 <- expression(p*3.141593*(d/2)^2)
DF7 <- cbind(d, p) 
RES7 <- propagate(expr = EXPR7, data = DF7, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES7
RES7$data
RES7$prop
RES7$sim
library(errors)
e <- set_errors(RES7$sim[1], RES7$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)


### Priklad
EXPR <- expression(a^b*x)
a = c(5, 0.1)
b = c(10, 0.1)
x = c(1, 0.1)
DAT <- cbind(a, b, x)
(res <- propagate(EXPR, DAT))
res$data
res$prop
res$sim
summary(res)
plot(res)


### Vypocítejte koeficient viskozity roztoku glycerinu Stokesovou metodou
r = c(0.0112, 0.0001) # polomer kulicky [cm]
l = c(31.23, 0.05) # draha kulicky za cas t [cm]
t = c(62.1, 0.2) # cas [s]
g = c(980.1,0) # tihove zrychleni [cm/s2] z tabulek
d0 = c(13.55,0) # hustota kulicky [g/cm3] z tabulek
d = c(1.28,0) # hustota roztoku [g/cm3] z tabulek
EXPR3 <- expression((2/9)*g*(((d0-d)*r^2)/l)*t)
DF3 <- cbind(r, l, t, g, d0, d) 
RES3 <- propagate(expr = EXPR3, data = DF3, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES3
RES3$data
RES3$prop
RES3$sim
EXPR4 <- expression((2/9)*980.1*(((13.55-1.28)*r^2)/l)*t)
DF4 <- cbind(r, l, t) 
RES4 <- propagate(expr = EXPR4, data = DF4, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES4$data
RES4$prop
RES4$sim
library(errors)
e <- set_errors(RES4$sim[1], RES4$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)


library(measurements) # prevod jednotek
### Vypocítejte koeficient viskozity roztoku glycerinu pomoci kapilarniho viskozimetru
p = c(20.12, 0.01) # tlak na vytoku z kapilary [mm Hg] to [Pa]
p = as.vector(conv_unit(c(20.12, 0.01), from="mmHg", to="Pa"))
r = c(0.570, 0.003) # polomer kapilary [mm] to [m]
r = conv_unit(c(0.570, 0.003), from="mm", to="m")
l = c(10.526, 0.005) # delka kapilary [mm] to [m]
l = conv_unit(c(10.526, 0.005), from="mm", to="m")
V = c(5.025, 0.001)# objem kapaliny vytekle za cas t [cm3] to [m3]
V = conv_unit(c(5.025, 0.001), from="cm3", to="m3")
t = c(27.34, 0.02) # cas [s]
pi
EXPR5 <- expression((3.141593*p*(r^4)*t)/(8*V*l))
DF5 <- cbind(p, r, l, V, t) 
RES5 <- propagate(expr = EXPR5, data = DF5, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES5 # [kg / m.s]
RES5$data
RES5$prop
RES5$sim
# conv_multiunit(x = RES5$sim[1:2], from="kg / m", to="g / cm")
library(errors)
e <- set_errors(RES5$sim[1], RES5$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)


## NEFUNGUJE pro jednu promennou, viz "metRology"
### Soucin rozpustnosti stríbrné soli AgX má hodnotu Ks = (4.0 +- 0.4) x 10-8. Jaká je chyba vypoctené rovnovázné koncentrace stríbrných iontu ve vode?
EXPR <- expression(sqrt(Ks))
Ks = c(4.0, 0.4)
DAT <- data.frame(Ks)
res <- propagate(EXPR,DAT)


##### Metoda 2: Výpocet pomocí intervalu 

### Priklad
EXPR3 <- expression(C * sqrt((520 * H * P)/(M *(t + 460))))
H <- c(64, 65)
M <- c(16, 16.2)
P <- c(361, 365)
t <- c(165, 170)
C <- c(38.4, 38.5)
DAT3 <- makeDat(EXPR3)
interval(DAT3, EXPR3, seq = 2)


### Priklad
EXPR5 <- expression(x^2 - x + 1)
x <- c(-2, 1)
curve(x^2 - x + 1, -2, 1)
DAT5 <- makeDat(EXPR5)
interval(DAT5, EXPR5, seq = 2)


library(metRology)

data(GUM.H.1)
GUM.H.1

## a simple uncertainty analysis for the product of two quantities
GUM(c("x1","x2"),c(2.3,1.1),c(0.030,0.015),c(5,9999),"x1*x2")

## a simple uncertainty analysis for the product of two quantities
GUM.validate(c("x1","x2"), c(2.3,1.1), c(0.030,0.015), c(5,9999), c("A","B"),c("Normal","Rectangular"),"x1*x2")


# Soucin rozpustnosti stríbrné soli AgX má hodnotu Ks = (4.0 +- 0.4) x 10-8. Jaká je chyba vypoctené rovnovázné koncentrace stríbrných iontu ve vode?
GUM("Ks",4.0e-8,0.4e-8,1,"sqrt(Ks)",sig.digits.U = 2) # tohle funguje
GUM.validate("Ks",4.0e-8,0.4e-8,1,"A","Normal","sqrt(Ks)")
GUM.validate("Ks",4.0e-8,0.4e-8,1,"B","Normal","sqrt(Ks)")

# Priklad
expr <- expression(a+b*2+c*3+d/2)
x <- list(a=1, b=3, c=2, d=11)
u <- lapply(x, function(x) x/10) # nejistota 10 %
u.expr<-uncert(expr, x, u, method="NUM")
u.expr
# function method
f <- function(a,b,c,d) a+b*2+c*3+d/2
u.fun<-uncert(f, x, u, method="NUM")
u.fun
# formula method
u.form<-uncert(~a+b*2+c*3+d/2, x, u, method="NUM")
u.form

# Priklad (s korelaci)
u.cor<-diag(1,4)
u.cor[3,4]<-u.cor[4,3]<-0.5
u.cor
# num
u.formc<-uncert(~a+b*2+c*3+d/2, x, u, method="NUM", cor=u.cor)
u.formc
# Monte Carlo
u.formc.MC<-uncert(~a+b*2+c*3+d/2, x, u, method="MC", cor=u.cor, B=200)
u.formc.MC


# Priklad
expr <- expression(a/(b-c))
x <- list(a=1, b=3, c=2)
u <- lapply(x, function(x) x/20)
set.seed(403)
u.invexpr<-uncertMC(expr, x, u, distrib=rep("norm", 3), B=999, keep.x=TRUE )
u.invexpr


library(errors)

x = c(0.9719378, 1.9840006, 2.9961830, 4.0123346, 5.0012799)
errors(x) = c(0.01, 0.01, 0.01, 0.01, 0.01)
y = c(0.9748992, 1.9627805, 2.9935831, 3.9921237, 4.9612555)
errors(y) = c(0.02, 0.02, 0.02, 0.02, 0.02)

# bez korelace
z <- x / y; z
#(z <- x / y)

# s korelací
correl(x, y) =  c(0.8864282, 0.4761841, -0.4140209, -0.4496266, 0.2911837) # zavedení korekace chyb
covar(x, y)
(z_correl <- x / y)