Category Archives: Allgemein

Basic R commands – 1.2 – permutations and matrix functions

05/05/2020Allgemein, Code, Programming, RMartin Stoppacher

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Basic Commands and Statistics with R - 1.2 - permutations
# path: ~/ownCloud/STA_Statistics/Basic_Statistics/
# file_name: statistic_basics3.R
# files_used:

#install.packages( 'gtools' )
library( gtools )
#install.packages( "combinat" )
library( combinat )

#install.packages( "magrittr" )
library( "magrittr" )

# - - - - - - - - - - - - - - - - 
#

x <- 6

factorial( x )
# [1] 720

for( i in 1:x ) print( factorial( i ) )
# [1] 1
# [1] 2
# [1] 6
# [1] 24
# [1] 120
# [1] 720

# - - - - - - - - - - - - - - - - 
#

x <- c( "a", "b", "c" )

combinat::permn( x )

# [[1]]
# [1] "a" "b" "c"
# [[2]]
# [1] "a" "c" "b"
# [[3]]
# [1] "c" "a" "b"
# [[4]]
# [1] "c" "b" "a"
# [[5]]
# [1] "b" "c" "a"
# [[6]]
# [1] "b" "a" "c"

# - - - - - - - - - - - - - - - - 
# string conversion using combinat::

x <- "ABCDEFG"
x <- strsplit( x, "" )[[1]]
y <- combinat::permn( x )
y

x <- LETTERS[1:5]
y <- combinat::permn( x )

paste( y[[1]], collapse = "" )
for( i in 1:length( y ) ) print( paste( y[[i]], collapse = "" ) )

# - - - - - - - - - - - - - - - - 
# using gtools:: for permutations

gtools::permutations( 5, 5, x )
# [,1] [,2] [,3] [,4] [,5]
# [1,] "A"  "C"  "I"  "M"  "P" 
# [2,] "A"  "C"  "I"  "P"  "M" 
# [3,] "A"  "C"  "M"  "I"  "P" 
# [4,] "A"  "C"  "M"  "P"  "I" 
# [5,] "A"  "C"  "P"  "I"  "M"
# ...

gtools::permutations( 5, 5, x, repeats.allowed = TRUE )
# [,1] [,2] [,3] [,4] [,5]
# [1,] "A"  "A"  "A"  "A"  "A" 
# [2,] "A"  "A"  "A"  "A"  "C" 
# [3,] "A"  "A"  "A"  "A"  "I" 
# [4,] "A"  "A"  "A"  "A"  "M" 
# [5,] "A"  "A"  "A"  "A"  "P" 
# [6,] "A"  "A"  "A"  "C"  "A" 
# [7,] "A"  "A"  "A"  "C"  "C" 
# [8,] "A"  "A"  "A"  "C"  "I" 
# [9,] "A"  "A"  "A"  "C"  "M" 
# [10,] "A"  "A"  "A"  "C"  "P" 
# ...

# - - - - - - - - - - - - - - - - 
# using gtools:: for combinations

gtools::combinations( 4, 2, x ) 

gtools::combinations( 5, 2, x )

gtools::combinations( 5, 3, x )

gtools::combinations( 4, 2, x, repeats.allowed = TRUE )

gtools::combinations( 2, 5, x, repeats.allowed = TRUE )

# Martin Stoppacher                                                             #
# office@martinstoppacher.com                                                   #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#################################################################################

100

101

102

103

104

105

106

107

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Basic Commands and Statistics with R - 1.2 - permutations

# path: ~/ownCloud/STA_Statistics/Basic_Statistics/

# file_name: statistic_basics3.R

# files_used:

#install.packages( 'gtools' )

library( gtools )

#install.packages( "combinat" )

library( combinat )

#install.packages( "magrittr" )

library( "magrittr" )

# - - - - - - - - - - - - - - - -

x <- 6

factorial( x )

# [1] 720

for( i in 1:x ) print( factorial( i ) )

# [1] 1

# [1] 2

# [1] 6

# [1] 24

# [1] 120

# [1] 720

# - - - - - - - - - - - - - - - -

x <- c( "a", "b", "c" )

combinat::permn( x )

# [[1]]

# [1] "a" "b" "c"

# [[2]]

# [1] "a" "c" "b"

# [[3]]

# [1] "c" "a" "b"

# [[4]]

# [1] "c" "b" "a"

# [[5]]

# [1] "b" "c" "a"

# [[6]]

# [1] "b" "a" "c"

# - - - - - - - - - - - - - - - -

# string conversion using combinat::

x <- "ABCDEFG"

x <- strsplit( x, "" )[[1]]

y <- combinat::permn( x )

x <- LETTERS[1:5]

y <- combinat::permn( x )

paste( y[[1]], collapse = "" )

for( i in 1:length( y ) ) print( paste( y[[i]], collapse = "" ) )

# - - - - - - - - - - - - - - - -

# using gtools:: for permutations

gtools::permutations( 5, 5, x )

# [,1] [,2] [,3] [,4] [,5]

# [1,] "A" "C" "I" "M" "P"

# [2,] "A" "C" "I" "P" "M"

# [3,] "A" "C" "M" "I" "P"

# [4,] "A" "C" "M" "P" "I"

# [5,] "A" "C" "P" "I" "M"

# ...

gtools::permutations( 5, 5, x, repeats.allowed = TRUE )

# [,1] [,2] [,3] [,4] [,5]

# [1,] "A" "A" "A" "A" "A"

# [2,] "A" "A" "A" "A" "C"

# [3,] "A" "A" "A" "A" "I"

# [4,] "A" "A" "A" "A" "M"

# [5,] "A" "A" "A" "A" "P"

# [6,] "A" "A" "A" "C" "A"

# [7,] "A" "A" "A" "C" "C"

# [8,] "A" "A" "A" "C" "I"

# [9,] "A" "A" "A" "C" "M"

# [10,] "A" "A" "A" "C" "P"

# ...

# - - - - - - - - - - - - - - - -

# using gtools:: for combinations

gtools::combinations( 4, 2, x )

gtools::combinations( 5, 2, x )

gtools::combinations( 5, 3, x )

gtools::combinations( 4, 2, x, repeats.allowed = TRUE )

gtools::combinations( 2, 5, x, repeats.allowed = TRUE )

# Martin Stoppacher #

# office@martinstoppacher.com #

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

#################################################################################

# matrix algebra

m1 <- matrix( c( 1, 2, 3, 4, 5, 6, 7, 8, 9 ) , 3, 3 )
m1

m2 <- matrix( rep( 2, 9 ), 3, 3 )
m2

#       [,1] [,2] [,3]
# [1,]    1    4    7
# [2,]    2    5    8
# [3,]    3    6    9

m1 * m1  # element multipliction

#       [,1] [,2] [,3]
# [1,]    1   16   49
# [2,]    4   25   64
# [3,]    9   36   81

m1 * m2

#       [,1] [,2] [,3]
# [1,]    2    8   14
# [2,]    4   10   16
# [3,]    6   12   18


# Matrix Multiplication
m1 %*% m1

#       [,1] [,2] [,3]
# [1,]   30   66  102
# [2,]   36   81  126
# [3,]   42   96  150

# Transpose
t( m1 )

#       [,1] [,2] [,3]
# [1,]    1    2    3
# [2,]    4    5    6
# [3,]    7    8    9

# Outer Product

m1[ 1 , ]
m1[ , 1 ] 
m1[ 1 , ] %o%  m1[ , 1 ] 

m1 %o% m1

# , , 1, 1
# 
#       [,1] [,2] [,3]
# [1,]    1    4    7
# [2,]    2    5    8
# [3,]    3    6    9
# 
# , , 2, 1
# 
#       [,1] [,2] [,3]
# [1,]    2    8   14
# [2,]    4   10   16
# [3,]    6   12   18
# 
# , , 3, 1
# 
#       [,1] [,2] [,3]
# [1,]    3   12   21
# [2,]    6   15   24
# [3,]    9   18   27
# 
# , , 1, 2
# 
#       [,1] [,2] [,3]
# [1,]    4   16   28
# [2,]    8   20   32
# [3,]   12   24   36
# 
# , , 2, 2
# 
#       [,1] [,2] [,3]
# [1,]    5   20   35
# [2,]   10   25   40
# [3,]   15   30   45
# 
# , , 3, 2
# 
#       [,1] [,2] [,3]
# [1,]    6   24   42
# [2,]   12   30   48
# [3,]   18   36   54
# 
# , , 1, 3
# 
#       [,1] [,2] [,3]
# [1,]    7   28   49
# [2,]   14   35   56
# [3,]   21   42   63
# 
# , , 2, 3
# 
#       [,1] [,2] [,3]
# [1,]    8   32   56
# [2,]   16   40   64
# [3,]   24   48   72
# 
# , , 3, 3
# 
#       [,1] [,2] [,3]
# [1,]    9   36   63
# [2,]   18   45   72
# [3,]   27   54   81

# kroenecker product

m1 %x% m1

#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,]    1    4    7    4   16   28    7   28   49
# [2,]    2    5    8    8   20   32   14   35   56
# [3,]    3    6    9   12   24   36   21   42   63
# [4,]    2    8   14    5   20   35    8   32   56
# [5,]    4   10   16   10   25   40   16   40   64
# [6,]    6   12   18   15   30   45   24   48   72
# [7,]    3   12   21    6   24   42    9   36   63
# [8,]    6   15   24   12   30   48   18   45   72
# [9,]    9   18   27   18   36   54   27   54   81

m1 %x% m2

#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,]    2    2    2    8    8    8   14   14   14
# [2,]    2    2    2    8    8    8   14   14   14
# [3,]    2    2    2    8    8    8   14   14   14
# [4,]    4    4    4   10   10   10   16   16   16
# [5,]    4    4    4   10   10   10   16   16   16
# [6,]    4    4    4   10   10   10   16   16   16
# [7,]    6    6    6   12   12   12   18   18   18
# [8,]    6    6    6   12   12   12   18   18   18
# [9,]    6    6    6   12   12   12   18   18   18

# cross Product 

m1[ , 1 ]
m1[ , 2 ] 

crossprod( m1, m2[ , 1 ]  )

#       [,1]
# [1,]   32
# [2,]   77
# [3,]  122

crossprod( m1 )

#      [,1] [,2] [,3]
# [1,]   14   32   50
# [2,]   32   77  122
# [3,]   50  122  194

crossprod( m1, m2 )

#       [,1] [,2] [,3]
# [1,]   12   12   12
# [2,]   30   30   30
# [3,]   48   48   48

# Diagonal Matrix with elements of x in the diagonal

diag( nrow = 3 )

# [,1] [,2] [,3]
# [1,]    1    0    0
# [2,]    0    1    0
# [3,]    0    0    1

diag( 3,  nrow = 3 )

# [,1] [,2] [,3]
# [1,]    3    0    0
# [2,]    0    3    0
# [3,]    0    0    3

# returns a vector of the principal diagonal elements
m1
diag( m1 )

# [,1] [,2] [,3]
# [1,]    1    0    0
# [2,]    0    2    0
# [3,]    0    0    3

# k * k identity matrix for m1[ , 1 ]
diag( m1[ , 1 ] )

# [,1] [,2] [,3]
# [1,]    1    0    0
# [2,]    0    2    0
# [3,]    0    0   

# simple combinations

cbind( m1, m2 )
rbind( m1, m2 )

# means and sums

rowMeans( m1 )
rowSums( m1 )

colMeans( m1 )
colSums( m1 )

# eigenvalues, eigenvectors

eigen( m1 )
e <- eigen( m1 )
e$values
e$vectors

# eigen() decomposition
# $values
# [1]  1.611684e+01 -1.116844e+00 -5.700691e-16
#
# $vectors
# [,1]       [,2]       [,3]
# [1,] -0.4645473 -0.8829060  0.4082483
# [2,] -0.5707955 -0.2395204 -0.8164966
# [3,] -0.6770438  0.4038651  0.4082483

# single value decomposition

svd( m1 )

# $d    singular values of m1
# [1] 1.684810e+01 1.068370e+00 5.543107e-16
#
# $u    left singular vectors of m1
# [,1]        [,2]       [,3]
# [1,] -0.4796712  0.77669099  0.4082483
# [2,] -0.5723678  0.07568647 -0.8164966
# [3,] -0.6650644 -0.62531805  0.4082483
#
# $v    right singular vectors of m1
# [,1]       [,2]       [,3]
# [1,] -0.2148372 -0.8872307  0.4082483
# [2,] -0.5205874 -0.2496440 -0.8164966
# [3,] -0.8263375  0.3879428  0.4082483


# QR decomposition 

qr( m1 )

# $qr
# [,1]      [,2]          [,3]
# [1,] -3.7416574 -8.552360 -1.336306e+01
# [2,]  0.5345225  1.963961  3.927922e+00
# [3,]  0.8017837  0.988693  1.776357e-15
#
# $rank
# [1] 2
#
# $qraux  
# [1] 1.267261e+00 1.149954e+00 1.776357e-15
#
# $pivot
# [1] 1 2 3
#
# attr(,"class")
# [1] "qr"


# inverse 

solve( m1 )

# Error in solve.default(m1) : 
#   Lapack routine dgesv: system is exactly singular: U[3,3] = 0

m3 <- matrix( c( 2, 6, 1, 4 ) , 2, 2 )
m3

# [,1] [,2]
# [1,]    2    1
# [2,]    6    4

sm3 <- solve( m3 )

# [,1] [,2]
# [1,]    2 -0.5
# [2,]   -3  1.0

m3 %*% sm3

# [,1] [,2]
# [1,]    1    0
# [2,]    0    1

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

# matrix algebra

m1 <- matrix( c( 1, 2, 3, 4, 5, 6, 7, 8, 9 ) , 3, 3 )

m2 <- matrix( rep( 2, 9 ), 3, 3 )

# [,1] [,2] [,3]

# [1,] 1 4 7

# [2,] 2 5 8

# [3,] 3 6 9

m1 * m1 # element multipliction

# [,1] [,2] [,3]

# [1,] 1 16 49

# [2,] 4 25 64

# [3,] 9 36 81

m1 * m2

# [,1] [,2] [,3]

# [1,] 2 8 14

# [2,] 4 10 16

# [3,] 6 12 18

# Matrix Multiplication

m1 %*% m1

# [,1] [,2] [,3]

# [1,] 30 66 102

# [2,] 36 81 126

# [3,] 42 96 150

# Transpose

t( m1 )

# [,1] [,2] [,3]

# [1,] 1 2 3

# [2,] 4 5 6

# [3,] 7 8 9

# Outer Product

m1[ 1 , ]

m1[ , 1 ]

m1[ 1 , ] %o% m1[ , 1 ]

m1 %o% m1

# , , 1, 1

# [,1] [,2] [,3]

# [1,] 1 4 7

# [2,] 2 5 8

# [3,] 3 6 9

# , , 2, 1

# [,1] [,2] [,3]

# [1,] 2 8 14

# [2,] 4 10 16

# [3,] 6 12 18

# , , 3, 1

# [,1] [,2] [,3]

# [1,] 3 12 21

# [2,] 6 15 24

# [3,] 9 18 27

# , , 1, 2

# [,1] [,2] [,3]

# [1,] 4 16 28

# [2,] 8 20 32

# [3,] 12 24 36

# , , 2, 2

# [,1] [,2] [,3]

# [1,] 5 20 35

# [2,] 10 25 40

# [3,] 15 30 45

# , , 3, 2

# [,1] [,2] [,3]

# [1,] 6 24 42

# [2,] 12 30 48

# [3,] 18 36 54

# , , 1, 3

# [,1] [,2] [,3]

# [1,] 7 28 49

# [2,] 14 35 56

# [3,] 21 42 63

# , , 2, 3

# [,1] [,2] [,3]

# [1,] 8 32 56

# [2,] 16 40 64

# [3,] 24 48 72

# , , 3, 3

# [,1] [,2] [,3]

# [1,] 9 36 63

# [2,] 18 45 72

# [3,] 27 54 81

# kroenecker product

m1 %x% m1

# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]

# [1,] 1 4 7 4 16 28 7 28 49

# [2,] 2 5 8 8 20 32 14 35 56

# [3,] 3 6 9 12 24 36 21 42 63

# [4,] 2 8 14 5 20 35 8 32 56

# [5,] 4 10 16 10 25 40 16 40 64

# [6,] 6 12 18 15 30 45 24 48 72

# [7,] 3 12 21 6 24 42 9 36 63

# [8,] 6 15 24 12 30 48 18 45 72

# [9,] 9 18 27 18 36 54 27 54 81

m1 %x% m2

# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]

# [1,] 2 2 2 8 8 8 14 14 14

# [2,] 2 2 2 8 8 8 14 14 14

# [3,] 2 2 2 8 8 8 14 14 14

# [4,] 4 4 4 10 10 10 16 16 16

# [5,] 4 4 4 10 10 10 16 16 16

# [6,] 4 4 4 10 10 10 16 16 16

# [7,] 6 6 6 12 12 12 18 18 18

# [8,] 6 6 6 12 12 12 18 18 18

# [9,] 6 6 6 12 12 12 18 18 18

# cross Product

m1[ , 1 ]

m1[ , 2 ]

crossprod( m1, m2[ , 1 ] )

# [,1]

# [1,] 32

# [2,] 77

# [3,] 122

crossprod( m1 )

# [,1] [,2] [,3]

# [1,] 14 32 50

# [2,] 32 77 122

# [3,] 50 122 194

crossprod( m1, m2 )

# [,1] [,2] [,3]

# [1,] 12 12 12

# [2,] 30 30 30

# [3,] 48 48 48

# Diagonal Matrix with elements of x in the diagonal

diag( nrow = 3 )

# [,1] [,2] [,3]

# [1,] 1 0 0

# [2,] 0 1 0

# [3,] 0 0 1

diag( 3, nrow = 3 )

# [,1] [,2] [,3]

# [1,] 3 0 0

# [2,] 0 3 0

# [3,] 0 0 3

# returns a vector of the principal diagonal elements

diag( m1 )

# [,1] [,2] [,3]

# [1,] 1 0 0

# [2,] 0 2 0

# [3,] 0 0 3

# k * k identity matrix for m1[ , 1 ]

diag( m1[ , 1 ] )

# [,1] [,2] [,3]

# [1,] 1 0 0

# [2,] 0 2 0

# [3,] 0 0

# simple combinations

cbind( m1, m2 )

rbind( m1, m2 )

# means and sums

rowMeans( m1 )

rowSums( m1 )

colMeans( m1 )

colSums( m1 )

# eigenvalues, eigenvectors

eigen( m1 )

e <- eigen( m1 )

e$values

e$vectors

# eigen() decomposition

# $values

# [1] 1.611684e+01 -1.116844e+00 -5.700691e-16

# $vectors

# [,1] [,2] [,3]

# [1,] -0.4645473 -0.8829060 0.4082483

# [2,] -0.5707955 -0.2395204 -0.8164966

# [3,] -0.6770438 0.4038651 0.4082483

# single value decomposition

svd( m1 )

# $d singular values of m1

# [1] 1.684810e+01 1.068370e+00 5.543107e-16

# $u left singular vectors of m1

# [,1] [,2] [,3]

# [1,] -0.4796712 0.77669099 0.4082483

# [2,] -0.5723678 0.07568647 -0.8164966

# [3,] -0.6650644 -0.62531805 0.4082483

# $v right singular vectors of m1

# [,1] [,2] [,3]

# [1,] -0.2148372 -0.8872307 0.4082483

# [2,] -0.5205874 -0.2496440 -0.8164966

# [3,] -0.8263375 0.3879428 0.4082483

# QR decomposition

qr( m1 )

# $qr

# [,1] [,2] [,3]

# [1,] -3.7416574 -8.552360 -1.336306e+01

# [2,] 0.5345225 1.963961 3.927922e+00

# [3,] 0.8017837 0.988693 1.776357e-15

# $rank

# [1] 2

# $qraux

# [1] 1.267261e+00 1.149954e+00 1.776357e-15

# $pivot

# [1] 1 2 3

# attr(,"class")

# [1] "qr"

# inverse

solve( m1 )

# Error in solve.default(m1) :

# Lapack routine dgesv: system is exactly singular: U[3,3] = 0

m3 <- matrix( c( 2, 6, 1, 4 ) , 2, 2 )

# [,1] [,2]

# [1,] 2 1

# [2,] 6 4

sm3 <- solve( m3 )

# [,1] [,2]

# [1,] 2 -0.5

# [2,] -3 1.0

m3 %*% sm3

# [,1] [,2]

# [1,] 1 0

# [2,] 0 1

Basic R commands – 1

04/21/2020Allgemein, Code, Data Analysis, Programming, R, StatisticsMartin Stoppacher

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Basic Commands and Statistics with R

rm(list = ls(all = TRUE))
getwd()
#system("ls")
setwd("~/ownCloud/STA_Statistics/basicR/")

search()

options(scipen=100)   # scientific off
options(scipen=0)     # scientific on
options(digits = 3)     

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# https://stat.ethz.ch/R-manual/R-patched/library/datasets/html/mtcars.html

attach(mtcars)
search()
detach(mtcars)
search()
mtcars         # The data was extracted from the 1974 Motor Trend US magazine 
# ?mtcars

class(mtcars)  # determine the class of an object

str(mtcars)    # Compactly display the internal structure of an R object

# 'data.frame':	32 obs. of  11 variables:

#                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
# Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
# Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
# Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
# Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
# Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
# Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
# Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
# Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
# Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
# Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
# Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
# Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
# Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
# Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
# Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
# Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
# Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
# Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
# Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
# Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
# Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
# Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
# AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
# Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
# Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
# Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
# Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
# Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
# Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
# Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
# Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
# Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

# A data frame with 32 observations on 11 (numeric) variables.

#[, 1]	mpg	Miles/(US) gallon
#[, 2]	cyl	Number of cylinders
#[, 3]	disp	Displacement (cu.in.)
#[, 4]	hp	Gross horsepower
#[, 5]	drat	Rear axle ratio
#[, 6]	wt	Weight (1000 lbs)
#[, 7]	qsec	1/4 mile time
#[, 8]	vs	Engine (0 = V-shaped, 1 = straight)
#[, 9]	am	Transmission (0 = automatic, 1 = manual)
#[,10]	gear	Number of forward gears
#[,11]	carb	Number of carburetors

colnames(mtcars) <- c("mpg","cyl","disp","hp","drat","wt",
                      "qsec","vs","am","gear","carb")

# colnames(mtcars) <- c("miles-per-gallon","cylinders","displacement",
#                       "horsepower","rear-axle-ratio","weight","qsec",
#                       "engine-vs","Transmissionam","gear","carburetors")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# data structures

letters
l <- letters
str(l)
# chr [1:26]  "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" 
#             "q" "r" "s" "t" "u" "v" "w" "x" "y" "z"

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # 
# Vectors

t <- c(1:3,"Hello",NA,FALSE,TRUE)
t
# [1] "1"     "2"     "3"     "Hello" NA      "FALSE" "TRUE" 
str(t)
# chr [1:7] "1" "2" "3" "Hello" NA "FALSE" "TRUE"

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Factors

mtcars$am[mtcars$am==1]
which(mtcars$am==1)

am.names <- mtcars$am
am.names[which(am.names==1)] <- "Automatic"
am.names
am.names[which(am.names==0)] <- "Manual"
am.names

am.names.factor <- factor(am.names)
am.names.factor
# [1] Automatic Automatic Automatic Manual    Manual  
# Levels: Automatic Manual

# internally stored as a table of:
# 1 Automatic
# 2 Manual

# Automatic and Manual are the levels of the factor 
levels(am.names.factor)

relevel(am.names.factor, ref="Manual")
# [1] Automatic Automatic Automatic Manual    Manual  
# Levels: Manual Automatic

am.names.factor[1]
as.character(am.names.factor[1])
as.numeric(am.names.factor[1])  # numeric takes the faktor not the level
levels(am.names.factor)[1]      # can be used to get the factor

table(am.names.factor)

is.factor(am.names.factor)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # 
# Data Frames

mtcars
str(mtcars)
# 'data.frame':	32 obs. of  11 variables:

df <- data.frame(1:20,60:41,letters[1:20])
df
colnames(df) <- c("n1","n2","letters")
df
t(df)

df$letters      # get the vectors from the dataframe
df$n1
df[,2]
df[1,3]
df[[3]]
levels(df$letters)
is.factor(df[,3])

df <- data.frame(1:20,60:41,letters[1:20], stringsAsFactors = FALSE)
df
is.factor(df[,3])

trees
mtcars
USMortality 
data()    # list of data associated with all current packages in the serch path

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Matrix

as.matrix(mtcars)
mtcars.mat <- as.matrix(mtcars) # all elements of a matrix have the same mode 
                                #(numeric, character)
x <- as.vector(mtcars.mat)      # all in order in the vector X
x

mat <- matrix(1:100,nco=10)
mat
as.vector(mat)

mat44 <- matrix(1:(4*4),nco=4)
mat44
as.vector(mat44)

dimnames(mtcars.mat)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Arrays

x <- 1:10000             # this is a matrix
dim(x) <- c(100,100)
x

x <- 1:1000              # array - a matrix with more than 2 dim
dim(x) <- c(10,10,10)
x

x <- 1:(2*3*3)           # array - a matrix with more than 2 dim
dim(x) <- c(2,3,3)
x

str(x)
# int [1:2, 1:3, 1:3] 1 2 3 4 5 6 7 8 9 10 ...

x <- letters[1:(2*3*3)]           # array - a matrix with more than 2 dim
dim(x) <- c(2,3,3)
x
str(x)
# chr [1:2, 1:3, 1:3] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" 
#                     "m" "n" "o" "p" "q" "r"

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Lists - collect different types of data objects

data()

l <- list(mtcars,barley,environmental,ethanol,melanoma)
                                            # most data from lattice package
str(l)

# List of 5
# $ :'data.frame':	32 obs. of  11 variables:
# ..$ mpg : num [1:32] 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
# ..$ cyl : num [1:32] 6 6 4 6 8 6 8 4 4 6 ...
# ..$ disp: num [1:32] 160 160 108 258 360 ...
# ..$ hp  : num [1:32] 110 110 93 110 175 105 245 62 95 123 ...
# ..$ drat: num [1:32] 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
# ..$ wt  : num [1:32] 2.62 2.88 2.32 3.21 3.44 ...
# ..$ qsec: num [1:32] 16.5 17 18.6 19.4 17 ...
# ..$ vs  : num [1:32] 0 0 1 1 0 1 0 1 1 1 ...
# ..$ am  : num [1:32] 1 1 1 0 0 0 0 0 0 0 ...
# ..$ gear: num [1:32] 4 4 4 3 3 3 3 4 4 4 ...
# ..$ carb: num [1:32] 4 4 1 1 2 1 4 2 2 4 ...
# $ :'data.frame':	120 obs. of  4 variables:
# ..$ yield  : num [1:120] 27 48.9 27.4 39.9 33 ...
# ..$ variety: Factor w/10 levels "Svansota","No.462",..: 3 3 3 3 3 3 7 7 7 7...
# ..$ year   : Factor w/2 levels "1932","1931": 2 2 2 2 2 2 2 2 2 2...
# ..$ site   : Factor w/6 levels "Grand Rapids",..: 3 6 4 5 1 2 3 6 4 5...
# $ :'data.frame':	111 obs. of  4 variables:
# ..$ ozone      : num [1:111] 41 36 12 18 23 19 8 16 11 14 ...
# ..$ radiation  : num [1:111] 190 118 149 313 299 99 19 256 290 274 ...
# ..$ temperature: num [1:111] 67 72 74 62 65 59 61 69 66 68 ...
# ..$ wind       : num [1:111] 7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...
# $ :'data.frame':	88 obs. of  3 variables:
# ..$ NOx: num [1:88] 3.74 2.29 1.5 2.88 0.76 ...
# ..$ C  : num [1:88] 12 12 12 12 12 9 9 9 12 12 ...
# ..$ E  : num [1:88] 0.907 0.761 1.108 1.016 1.189 ...
# $ :'data.frame':	37 obs. of  2 variables:
# ..$ year     : num [1:37] 1936 1937 1938 1939 1940 ...
# ..$ incidence: num [1:37] 0.9 0.8 0.8 1.3 1.4 1.2 1.7 1.8 1.6 1.5 ...

l[[1]]
l[[2]]
l[[3]]

l[[1]]

model <- lm(mpg~wt,data=l[[1]])   # building a model from a list object
names(model) 
str(model)                        # the model itselfe is returned as a list

model$coefficients
model[[1]]

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# some basic R functions

mtcars$mpg
# [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 ...
mtcars$mpg^2
# [1]  441.00  441.00 519.84 457.96 349.69 327.61 204.49  595.36  519.84  ...
sqrt(mtcars$mpg^2)
# [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 ...

pi
# [1] 3.141593
90*pi/180 # radiant conversion
sin(90*pi/180)
cos(90*pi/180)
cos(0*pi/180)

celsius <- 20
9/5*celsius+32

celsius <- -20:40
fahrenheit <- 9/5*celsius+32
plot(celsius,fahrenheit,type="s")

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Basic Commands and Statistics with R

rm(list = ls(all = TRUE))

getwd()

#system("ls")

setwd("~/ownCloud/STA_Statistics/basicR/")

search()

options(scipen=100) # scientific off

options(scipen=0) # scientific on

options(digits = 3)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# https://stat.ethz.ch/R-manual/R-patched/library/datasets/html/mtcars.html

attach(mtcars)

search()

detach(mtcars)

search()

mtcars # The data was extracted from the 1974 Motor Trend US magazine

# ?mtcars

class(mtcars) # determine the class of an object

str(mtcars) # Compactly display the internal structure of an R object

# 'data.frame': 32 obs. of 11 variables:

# mpg cyl disp hp drat wt qsec vs am gear carb

# Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4

# Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4

# Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1

# Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1

# Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2

# Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1

# Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4

# Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2

# Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2

# Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4

# Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4

# Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3

# Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3

# Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3

# Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4

# Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4

# Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4

# Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1

# Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2

# Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1

# Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1

# Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2

# AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2

# Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4

# Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2

# Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1

# Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2

# Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2

# Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4

# Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6

# Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8

# Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2

# A data frame with 32 observations on 11 (numeric) variables.

#[, 1] mpg Miles/(US) gallon

#[, 2] cyl Number of cylinders

#[, 3] disp Displacement (cu.in.)

#[, 4] hp Gross horsepower

#[, 5] drat Rear axle ratio

#[, 6] wt Weight (1000 lbs)

#[, 7] qsec 1/4 mile time

#[, 8] vs Engine (0 = V-shaped, 1 = straight)

#[, 9] am Transmission (0 = automatic, 1 = manual)

#[,10] gear Number of forward gears

#[,11] carb Number of carburetors

colnames(mtcars) <- c("mpg","cyl","disp","hp","drat","wt",

"qsec","vs","am","gear","carb")

# colnames(mtcars) <- c("miles-per-gallon","cylinders","displacement",

# "horsepower","rear-axle-ratio","weight","qsec",

# "engine-vs","Transmissionam","gear","carburetors")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# data structures

letters

l <- letters

str(l)

# chr [1:26] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p"

# "q" "r" "s" "t" "u" "v" "w" "x" "y" "z"

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Vectors

t <- c(1:3,"Hello",NA,FALSE,TRUE)

# [1] "1" "2" "3" "Hello" NA "FALSE" "TRUE"

str(t)

# chr [1:7] "1" "2" "3" "Hello" NA "FALSE" "TRUE"

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Factors

mtcars$am[mtcars$am==1]

which(mtcars$am==1)

am.names <- mtcars$am

am.names[which(am.names==1)] <- "Automatic"

am.names

am.names[which(am.names==0)] <- "Manual"

am.names

am.names.factor <- factor(am.names)

am.names.factor

# [1] Automatic Automatic Automatic Manual Manual

# Levels: Automatic Manual

# internally stored as a table of:

# 1 Automatic

# 2 Manual

# Automatic and Manual are the levels of the factor

levels(am.names.factor)

relevel(am.names.factor, ref="Manual")

# [1] Automatic Automatic Automatic Manual Manual

# Levels: Manual Automatic

am.names.factor[1]

as.character(am.names.factor[1])

as.numeric(am.names.factor[1]) # numeric takes the faktor not the level

levels(am.names.factor)[1] # can be used to get the factor

table(am.names.factor)

is.factor(am.names.factor)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Data Frames

mtcars

str(mtcars)

# 'data.frame': 32 obs. of 11 variables:

df <- data.frame(1:20,60:41,letters[1:20])

colnames(df) <- c("n1","n2","letters")

t(df)

df$letters # get the vectors from the dataframe

df$n1

df[,2]

df[1,3]

df[[3]]

levels(df$letters)

is.factor(df[,3])

df <- data.frame(1:20,60:41,letters[1:20], stringsAsFactors = FALSE)

is.factor(df[,3])

trees

mtcars

USMortality

data() # list of data associated with all current packages in the serch path

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Matrix

as.matrix(mtcars)

mtcars.mat <- as.matrix(mtcars) # all elements of a matrix have the same mode

#(numeric, character)

x <- as.vector(mtcars.mat) # all in order in the vector X

mat <- matrix(1:100,nco=10)

mat

as.vector(mat)

mat44 <- matrix(1:(4*4),nco=4)

mat44

as.vector(mat44)

dimnames(mtcars.mat)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Arrays

x <- 1:10000 # this is a matrix

dim(x) <- c(100,100)

x <- 1:1000 # array - a matrix with more than 2 dim

dim(x) <- c(10,10,10)

x <- 1:(2*3*3) # array - a matrix with more than 2 dim

dim(x) <- c(2,3,3)

str(x)

# int [1:2, 1:3, 1:3] 1 2 3 4 5 6 7 8 9 10 ...

x <- letters[1:(2*3*3)] # array - a matrix with more than 2 dim

dim(x) <- c(2,3,3)

str(x)

# chr [1:2, 1:3, 1:3] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l"

# "m" "n" "o" "p" "q" "r"

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Lists - collect different types of data objects

data()

l <- list(mtcars,barley,environmental,ethanol,melanoma)

# most data from lattice package

str(l)

# List of 5

# $ :'data.frame': 32 obs. of 11 variables:

# ..$ mpg : num [1:32] 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...

# ..$ cyl : num [1:32] 6 6 4 6 8 6 8 4 4 6 ...

# ..$ disp: num [1:32] 160 160 108 258 360 ...

# ..$ hp : num [1:32] 110 110 93 110 175 105 245 62 95 123 ...

# ..$ drat: num [1:32] 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...

# ..$ wt : num [1:32] 2.62 2.88 2.32 3.21 3.44 ...

# ..$ qsec: num [1:32] 16.5 17 18.6 19.4 17 ...

# ..$ vs : num [1:32] 0 0 1 1 0 1 0 1 1 1 ...

# ..$ am : num [1:32] 1 1 1 0 0 0 0 0 0 0 ...

# ..$ gear: num [1:32] 4 4 4 3 3 3 3 4 4 4 ...

# ..$ carb: num [1:32] 4 4 1 1 2 1 4 2 2 4 ...

# $ :'data.frame': 120 obs. of 4 variables:

# ..$ yield : num [1:120] 27 48.9 27.4 39.9 33 ...

# ..$ variety: Factor w/10 levels "Svansota","No.462",..: 3 3 3 3 3 3 7 7 7 7...

# ..$ year : Factor w/2 levels "1932","1931": 2 2 2 2 2 2 2 2 2 2...

# ..$ site : Factor w/6 levels "Grand Rapids",..: 3 6 4 5 1 2 3 6 4 5...

# $ :'data.frame': 111 obs. of 4 variables:

# ..$ ozone : num [1:111] 41 36 12 18 23 19 8 16 11 14 ...

# ..$ radiation : num [1:111] 190 118 149 313 299 99 19 256 290 274 ...

# ..$ temperature: num [1:111] 67 72 74 62 65 59 61 69 66 68 ...

# ..$ wind : num [1:111] 7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...

# $ :'data.frame': 88 obs. of 3 variables:

# ..$ NOx: num [1:88] 3.74 2.29 1.5 2.88 0.76 ...

# ..$ C : num [1:88] 12 12 12 12 12 9 9 9 12 12 ...

# ..$ E : num [1:88] 0.907 0.761 1.108 1.016 1.189 ...

# $ :'data.frame': 37 obs. of 2 variables:

# ..$ year : num [1:37] 1936 1937 1938 1939 1940 ...

# ..$ incidence: num [1:37] 0.9 0.8 0.8 1.3 1.4 1.2 1.7 1.8 1.6 1.5 ...

l[[1]]

l[[2]]

l[[3]]

l[[1]]

model <- lm(mpg~wt,data=l[[1]]) # building a model from a list object

names(model)

str(model) # the model itselfe is returned as a list

model$coefficients

model[[1]]

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# some basic R functions

mtcars$mpg

# [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 ...

mtcars$mpg^2

# [1] 441.00 441.00 519.84 457.96 349.69 327.61 204.49 595.36 519.84 ...

sqrt(mtcars$mpg^2)

# [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2 ...

# [1] 3.141593

90*pi/180 # radiant conversion

sin(90*pi/180)

cos(90*pi/180)

cos(0*pi/180)

celsius <- 20

9/5*celsius+32

celsius <- -20:40

fahrenheit <- 9/5*celsius+32

plot(celsius,fahrenheit,type="s")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
sum(mtcars$mpg)
# [1] 642.9
cumsum(mtcars$mpg)
# [1]  21.0  42.0  64.8  86.2 104.9 123.0 137.3 161.7
cumprod(mtcars$mpg)
# [1] 2.100000e+01 4.410000e+02 1.005480e+04 2.151727e+05 
#     4.023730e+06 7.282951e+07

c(1,2,3,4,5)    # vectors # concatenating function
c(1:5)
c(TRUE,FALSE)
c("TRUE","FALSE") # character vector

x <- runif(5)
x
# [1] 0.832 0.101 0.926 0.253 0.619
sort(x)

order(x)
# [1] 0.101 0.253 0.619 0.832 0.926
x[order(x)]
# [1] 0.101 0.253 0.619 0.832 0.926

x <- c(1:5)
x
y <- c("TRUE","FALSE",NA) 
y
z <- c(x,y)
z
# [1] "1"     "2"     "3"     "4"     "5"     "TRUE"  "FALSE" NA    
rev(z)
# [1] NA      "FALSE" "TRUE"  "5"     "4"     "3"     "2"     "1"   

z[2]
z[2:5]
# [1] "2" "3" "4" "5"
z[z<4]
# [1] "1" "2" "3" NA 
z[z>4]
# [1] "5"     "TRUE"  "FALSE" NA 
z[z>=4]
# [1] "4"     "5"     "TRUE"  "FALSE" NA  
z[-c(3,4)]
# [1] "1"     "2"     "5"     "TRUE"  "FALSE" NA 
z[c(3,4)]
# [1] "3" "4"

z[is.na(z)]
z[is.na(z)]<-0
z
z[z==TRUE]<-1
z[z==FALSE]<-0
z
# [1] "1" "2" "3" "4" "5" "1" "0" "0"

plot(z)
plot(z,type="b")

x <- runif(length(z))*5

plot(x,z,type="b")
plot(x~z,type="b")

plot(z,x,type="b")
lines(z~x,type="b",col="red")

plot(sin((1:360)*pi/180),type="l")
plot(sin((1:360)*pi/180),cos((1:360)*pi/180),type="l")

plot(sin((1:360)*pi/18),cos((1:360)*pi/10),type="l")

plot(sin((1:360)*pi/18),cos((1:360)*pi/10),type="l")

plot(tan((1:360)*pi/180),cos((1:360)*pi/180),type="l")
plot(tan((1:360)*pi/180),sin((1:360)*pi/180),type="l")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

sum(mtcars$mpg)

# [1] 642.9

cumsum(mtcars$mpg)

# [1] 21.0 42.0 64.8 86.2 104.9 123.0 137.3 161.7

cumprod(mtcars$mpg)

# [1] 2.100000e+01 4.410000e+02 1.005480e+04 2.151727e+05

# 4.023730e+06 7.282951e+07

c(1,2,3,4,5) # vectors # concatenating function

c(1:5)

c(TRUE,FALSE)

c("TRUE","FALSE") # character vector

x <- runif(5)

# [1] 0.832 0.101 0.926 0.253 0.619

sort(x)

order(x)

# [1] 0.101 0.253 0.619 0.832 0.926

x[order(x)]

# [1] 0.101 0.253 0.619 0.832 0.926

x <- c(1:5)

y <- c("TRUE","FALSE",NA)

z <- c(x,y)

# [1] "1" "2" "3" "4" "5" "TRUE" "FALSE" NA

rev(z)

# [1] NA "FALSE" "TRUE" "5" "4" "3" "2" "1"

z[2]

z[2:5]

# [1] "2" "3" "4" "5"

z[z<4]

# [1] "1" "2" "3" NA

z[z>4]

# [1] "5" "TRUE" "FALSE" NA

z[z>=4]

# [1] "4" "5" "TRUE" "FALSE" NA

z[-c(3,4)]

# [1] "1" "2" "5" "TRUE" "FALSE" NA

z[c(3,4)]

# [1] "3" "4"

z[is.na(z)]

z[is.na(z)]<-0

z[z==TRUE]<-1

z[z==FALSE]<-0

# [1] "1" "2" "3" "4" "5" "1" "0" "0"

plot(z)

plot(z,type="b")

x <- runif(length(z))*5

plot(x,z,type="b")

plot(x~z,type="b")

plot(z,x,type="b")

lines(z~x,type="b",col="red")

plot(sin((1:360)*pi/180),type="l")

plot(sin((1:360)*pi/180),cos((1:360)*pi/180),type="l")

plot(sin((1:360)*pi/18),cos((1:360)*pi/10),type="l")

plot(tan((1:360)*pi/180),cos((1:360)*pi/180),type="l")

plot(tan((1:360)*pi/180),sin((1:360)*pi/180),type="l")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# some loops

for(i in 1:360){
    plot(cos((1:360)*pi/i),sin((1:360)*pi/i),type="l")
    print(i)
    Sys.sleep(0.1)
}

for(i in 360:1){
  plot(cos((1:360)*pi/i),sin((1:360)*pi/i),type="l")
  print(i)
  Sys.sleep(0.1)
}

for(i in 1:360){
  for(e in 1:360){
  plot(cos((1:360)*pi/i),sin((1:360)*pi/e),type="l")
  print(i)
  Sys.sleep(0.1)
  }
}

for(i in 1:360){
  for(e in 1:360){
    plot(sin((1:360)*pi/i),cos((1:360)*pi/e),type="l")
    print(i)
    Sys.sleep(0.1)
  }
}

for(i in 1:10){
  for(e in 1:10){
    plot(sin((1:360)*pi/i),cos((1:360)*pi/e),type="l")
    mtext(paste("i:",i,"e:",e), side=3, outer=TRUE, line=-3)
    print(paste("i:",i,"e:",e))
    Sys.sleep(0.1)
  }
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# function collection

letters
# "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" 
# "s" "t" "u" "v" "w" "x" "y" "z"
LETTERS
# "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" 
# "S" "T" "U" "V" "W" "X" "Y" "Z"
letters[1:15]
# [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o"

plot(mtcars$mpg)
print(mtcars$mpg)

table(mtcars$mpg)
table(mtcars$cyl)
#   4   6   8 
#   11  7   14 

table(mtcars[,9:10])
str(table(mtcars[,9:10])) # table with 2 dimmensions
table(mtcars[,9:11])
str(table(mtcars[,9:11]))
# 'table' int [1:2, 1:3, 1:6] 3 0 0 4 0 0 4 0 2 2 ...
# - attr(*, "dimnames")=List of 3
# ..$ am  : chr [1:2] "0" "1"
# ..$ gear: chr [1:3] "3" "4" "5"
# ..$ carb: chr [1:6] "1" "2" "3" "4" ...

length(mtcars$mpg)
cat(mtcars$mpg)
mean(mtcars$mpg)
median(mtcars$mpg)
range(mtcars$mpg)
unique(mtcars$mpg)

rep(mtcars$mpg,10)

names()
colnames()
rownames()

diff(mtcars$mpg)
plot(diff(mtcars$mpg),type="h",xlab="",ylab="difference")
points(diff(mtcars$mpg),col="red")

sort(mtcars$mpg)
order(mtcars$mpg)
rev(mtcars$mpg)
rev(sort(mtcars$mpg))

cumsum(mtcars$mpg)
cumprod(mtcars$mpg)

rank(mtcars$vs) # Returns the sample ranks of the values in a vector.
(r1 <- rank(x1 <- c(3, 1, 4, 15, 92)))
rank(mtcars$vs, ties.method= "first")  # first occurrence wins 
                                       ## ranks without averaging
rank(mtcars$vs, ties.method= "last")   #  last occurrence wins 
                                       ## ranks without averaging
rank(mtcars$vs, ties.method= "random") # ties broken at random 
                                       ## ranks without averaging
rank(mtcars$vs, ties.method= "random") # and again 
                                       ## ranks without averaging

1:10
7:20
intersect(1:10, 7:20)

match(1:10,7:20)    # match returns a vector of the positions of (first) 
                    # matches of its first argument in its second.

1:10 %in% c(1,3,5,9)
sstr <- c("c","ab","B","bba","c",NA,"@","bla","a","Ba","%")
sstr %in% c(letters, LETTERS)
sstr[sstr %in% c(letters, LETTERS)]

apply()
sapply()
aggregate()
aggregate(state.x77, list(Region = state.region), mean)
tapply()

merge(mtcars$vs,c(1,2))

read.csv()
read.table()

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# histograms

hist(mtcars$mpg)

# looping with R
par(mfrow=c(3,4))
for(i in 1:length(mtcars[1,])){ 
  hist(mtcars[,i],main=paste("Data: MTcars",colnames(mtcars)[i]),
                              xlab=paste(colnames(mtcars)[i]))
}
par(mfrow=c(1,1))

hist(mtcars$mpg, main="Data: MTcars - mpg Miles per gallon", 
     xlab="mpg Miles per gallon")

hist(mtcars$mpg, breaks=10, main="Data: MTcars - mpg Miles per gallon", 
     xlab="mpg Miles per gallon")

hist(mtcars$mpg, breaks=length(mtcars$mpg), 
     main="Data: MTcars - mpg Miles per gallon", 
     xlab="mpg Miles per gallon")

# probability densities, component
hist(mtcars$mpg, breaks=length(mtcars$mpg), freq=FALSE, 
     main="Data: MTcars - mpg Miles per gallon", xlab="mpg Miles per gallon")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#  Boxplot with Scatterplot

par( fig = c( 0, 0.8, 0, 0.8 ), new = TRUE)
plot( mtcars$wt, mtcars$mpg, xlab = "Car Weight",
      ylab = "Miles Per Gallon" )
abline( lm( mpg~wt ), col = "red" )      # regression line (y~x)
lines( lowess( wt, mpg ), col = "blue" ) # lowess line (x,y) 
                                         # LOWESS smoother which uses 
                                         # locally-weighted 
                                         # polynomial regression
model <- lm( mpg~wt )
text( 4, 32, paste( "Intercept", model$coefficients[1] ) )
text( 4, 30, paste( "wt", model$coefficients[2] ) )

par( fig = c( 0, 0.8, 0.55, 1 ), new = TRUE )
boxplot(mtcars$wt, horizontal = TRUE, axes=FALSE, 
                                 col= "green", notch = TRUE )
m.wt <- mean( mtcars$wt )
mtext( paste( "mean wt", m.wt ), side = 3, outer = TRUE, line = -8 )

par( fig = c( 0.65, 1, 0, 0.8 ), new = TRUE)
boxplot( mtcars$mpg, axes = FALSE, col = "green", notch = TRUE)
m.mpg <- mean( mtcars$mpg )
mtext( paste( "mean mpg", m.mpg ), side = 2, outer = TRUE, line = -28)

mtext( "Miles per gallon vs Car Weight", side = 3, outer = TRUE, line = -3 )

par( mfrow = c( 1, 1 ) )

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# some loops

for(i in 1:360){

plot(cos((1:360)*pi/i),sin((1:360)*pi/i),type="l")

print(i)

Sys.sleep(0.1)

}

for(i in 360:1){

plot(cos((1:360)*pi/i),sin((1:360)*pi/i),type="l")

print(i)

Sys.sleep(0.1)

}

for(i in 1:360){

for(e in 1:360){

plot(cos((1:360)*pi/i),sin((1:360)*pi/e),type="l")

print(i)

Sys.sleep(0.1)

}

for(i in 1:360){

for(e in 1:360){

plot(sin((1:360)*pi/i),cos((1:360)*pi/e),type="l")

print(i)

Sys.sleep(0.1)

}

for(i in 1:10){

for(e in 1:10){

plot(sin((1:360)*pi/i),cos((1:360)*pi/e),type="l")

mtext(paste("i:",i,"e:",e), side=3, outer=TRUE, line=-3)

print(paste("i:",i,"e:",e))

Sys.sleep(0.1)

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# function collection

letters

# "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r"

# "s" "t" "u" "v" "w" "x" "y" "z"

LETTERS

# "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R"

# "S" "T" "U" "V" "W" "X" "Y" "Z"

letters[1:15]

# [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o"

plot(mtcars$mpg)

print(mtcars$mpg)

table(mtcars$mpg)

table(mtcars$cyl)

# 4 6 8

# 11 7 14

table(mtcars[,9:10])

str(table(mtcars[,9:10])) # table with 2 dimmensions

table(mtcars[,9:11])

str(table(mtcars[,9:11]))

# 'table' int [1:2, 1:3, 1:6] 3 0 0 4 0 0 4 0 2 2 ...

# - attr(*, "dimnames")=List of 3

# ..$ am : chr [1:2] "0" "1"

# ..$ gear: chr [1:3] "3" "4" "5"

# ..$ carb: chr [1:6] "1" "2" "3" "4" ...

length(mtcars$mpg)

cat(mtcars$mpg)

mean(mtcars$mpg)

median(mtcars$mpg)

range(mtcars$mpg)

unique(mtcars$mpg)

rep(mtcars$mpg,10)

names()

colnames()

rownames()

diff(mtcars$mpg)

plot(diff(mtcars$mpg),type="h",xlab="",ylab="difference")

points(diff(mtcars$mpg),col="red")

sort(mtcars$mpg)

order(mtcars$mpg)

rev(mtcars$mpg)

rev(sort(mtcars$mpg))

cumsum(mtcars$mpg)

cumprod(mtcars$mpg)

rank(mtcars$vs) # Returns the sample ranks of the values in a vector.

(r1 <- rank(x1 <- c(3, 1, 4, 15, 92)))

rank(mtcars$vs, ties.method= "first") # first occurrence wins

## ranks without averaging

rank(mtcars$vs, ties.method= "last") # last occurrence wins

## ranks without averaging

rank(mtcars$vs, ties.method= "random") # ties broken at random

## ranks without averaging

rank(mtcars$vs, ties.method= "random") # and again

## ranks without averaging

1:10

7:20

intersect(1:10, 7:20)

match(1:10,7:20) # match returns a vector of the positions of (first)

# matches of its first argument in its second.

1:10 %in% c(1,3,5,9)

sstr <- c("c","ab","B","bba","c",NA,"@","bla","a","Ba","%")

sstr %in% c(letters, LETTERS)

sstr[sstr %in% c(letters, LETTERS)]

apply()

sapply()

aggregate()

aggregate(state.x77, list(Region = state.region), mean)

tapply()

merge(mtcars$vs,c(1,2))

read.csv()

read.table()

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# histograms

hist(mtcars$mpg)

# looping with R

par(mfrow=c(3,4))

for(i in 1:length(mtcars[1,])){

hist(mtcars[,i],main=paste("Data: MTcars",colnames(mtcars)[i]),

xlab=paste(colnames(mtcars)[i]))

}

par(mfrow=c(1,1))

hist(mtcars$mpg, main="Data: MTcars - mpg Miles per gallon",

xlab="mpg Miles per gallon")

hist(mtcars$mpg, breaks=10, main="Data: MTcars - mpg Miles per gallon",

xlab="mpg Miles per gallon")

hist(mtcars$mpg, breaks=length(mtcars$mpg),

main="Data: MTcars - mpg Miles per gallon",

xlab="mpg Miles per gallon")

# probability densities, component

hist(mtcars$mpg, breaks=length(mtcars$mpg), freq=FALSE,

main="Data: MTcars - mpg Miles per gallon", xlab="mpg Miles per gallon")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Boxplot with Scatterplot

par( fig = c( 0, 0.8, 0, 0.8 ), new = TRUE)

plot( mtcars$wt, mtcars$mpg, xlab = "Car Weight",

ylab = "Miles Per Gallon" )

abline( lm( mpg~wt ), col = "red" ) # regression line (y~x)

lines( lowess( wt, mpg ), col = "blue" ) # lowess line (x,y)

# LOWESS smoother which uses

# locally-weighted

# polynomial regression

model <- lm( mpg~wt )

text( 4, 32, paste( "Intercept", model$coefficients[1] ) )

text( 4, 30, paste( "wt", model$coefficients[2] ) )

par( fig = c( 0, 0.8, 0.55, 1 ), new = TRUE )

boxplot(mtcars$wt, horizontal = TRUE, axes=FALSE,

col= "green", notch = TRUE )

m.wt <- mean( mtcars$wt )

mtext( paste( "mean wt", m.wt ), side = 3, outer = TRUE, line = -8 )

par( fig = c( 0.65, 1, 0, 0.8 ), new = TRUE)

boxplot( mtcars$mpg, axes = FALSE, col = "green", notch = TRUE)

m.mpg <- mean( mtcars$mpg )

mtext( paste( "mean mpg", m.mpg ), side = 2, outer = TRUE, line = -28)

mtext( "Miles per gallon vs Car Weight", side = 3, outer = TRUE, line = -3 )

par( mfrow = c( 1, 1 ) )

# Martin Stoppacher                                                             #
# office@martinstoppacher.com                                                   #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#################################################################################

# Martin Stoppacher #

# office@martinstoppacher.com #

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

#################################################################################

Renewable electricity output – Regions

12/07/2019Allgemein, Data Analysis, Programming, RMartin Stoppacher

https://host6.ssl-net.net/martinstoppacher_com/blog/?p=1194

# -------------

# region

# World
# Euro area
# European Union
# Europe & Central Asia
# Arab World
# North America
# Latin America & Caribbean
# Middle East & North Africa
# Sub-Saharan Africa 
# East Asia & Pacific
# South Asia

target.area <- which(renewable.electricity.output$country=="World")
renewable.electricity.output[target.area,]
plot(renewable.electricity.output[target.area,4:3],type="l",
     main=titel,sub=subtitel,ylab="%",xlab="",xlim=c(1980,2020),ylim=c(1,90),
     axes=TRUE)
text(1985, 20, "World",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Euro area")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="red")
text(1985, 12, "Euro area", col="red",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Arab World")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="green")
text(1985, 5, "Arab World", col="green",cex=0.7)

target.area <- which(renewable.electricity.output$country=="North America")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="blue",lty=2)
text(1985, 17, "North America", col="blue",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Latin America & Caribbean")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="yellow")
text(1985, 65, "Latin America & Caribbean", col="yellow",cex=0.7)

# target.area <- which(renewable.electricity.output$country=="Middle East & North Africa")
# renewable.electricity.output[target.area,]
# lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="orange")
# text(1985, 5, "Middle East & North Africa", col="orange",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Sub-Saharan Africa")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="cyan",lty=3)
text(1985, 23, "Sub-Saharan Africa", col="cyan",cex=0.7)

target.area <- which(renewable.electricity.output$country=="East Asia & Pacific")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="magenta",lty=4)
text(1985, 14.5, "East Asia & Pacific", col="magenta",cex=0.7)

target.area <- which(renewable.electricity.output$country=="South Asia")
renewable.electricity.output[target.area,]
lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="pink",lty=5)
text(1985, 28, "South Asia", col="pink",cex=0.7)

# -------------

# region

# World

# Euro area

# European Union

# Europe & Central Asia

# Arab World

# North America

# Latin America & Caribbean

# Middle East & North Africa

# Sub-Saharan Africa

# East Asia & Pacific

# South Asia

target.area <- which(renewable.electricity.output$country=="World")

renewable.electricity.output[target.area,]

plot(renewable.electricity.output[target.area,4:3],type="l",

main=titel,sub=subtitel,ylab="%",xlab="",xlim=c(1980,2020),ylim=c(1,90),

axes=TRUE)

text(1985, 20, "World",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Euro area")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="red")

text(1985, 12, "Euro area", col="red",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Arab World")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="green")

text(1985, 5, "Arab World", col="green",cex=0.7)

target.area <- which(renewable.electricity.output$country=="North America")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="blue",lty=2)

text(1985, 17, "North America", col="blue",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Latin America & Caribbean")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="yellow")

text(1985, 65, "Latin America & Caribbean", col="yellow",cex=0.7)

# target.area <- which(renewable.electricity.output$country=="Middle East & North Africa")

# renewable.electricity.output[target.area,]

# lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="orange")

# text(1985, 5, "Middle East & North Africa", col="orange",cex=0.7)

target.area <- which(renewable.electricity.output$country=="Sub-Saharan Africa")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="cyan",lty=3)

text(1985, 23, "Sub-Saharan Africa", col="cyan",cex=0.7)

target.area <- which(renewable.electricity.output$country=="East Asia & Pacific")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="magenta",lty=4)

text(1985, 14.5, "East Asia & Pacific", col="magenta",cex=0.7)

target.area <- which(renewable.electricity.output$country=="South Asia")

renewable.electricity.output[target.area,]

lines(renewable.electricity.output[target.area,4:3],type="l",main=titel,ylab="%",xlab="",col="pink",lty=5)

text(1985, 28, "South Asia", col="pink",cex=0.7)

Renewable electricity output – World Map

12/07/2019Allgemein, Data Analysis, Programming, RMartin Stoppacher

https://host6.ssl-net.net/martinstoppacher_com/blog/?p=1194

# ----------------

library(rworldmap)

options("scipen"=100, "digits"=4)

year <- 2015
target.year <- which(renewable.electricity.output$year==year)
head(renewable.electricity.output[target.year,])

renewable.electricity.output.year <- na.omit(renewable.electricity.output[target.year,])
renewable.electricity.output.year <- renewable.electricity.output.year[which(renewable.electricity.output.year$EG.ELC.RNEW.ZS>0),]

mapped_data <- joinCountryData2Map(renewable.electricity.output.year, joinCode = "ISO2", nameJoinColumn = "iso2c")
par(mai=c(0,0,0.2,0),xaxs="i",yaxs="i")
mapCountryData(mapped_data, nameColumnToPlot = "EG.ELC.RNEW.ZS",mapTitle=paste(titel,subtitel))

mapCountryData(mapped_data, nameColumnToPlot = "EG.ELC.RNEW.ZS",mapRegion="world",mapTitle=paste(titel,subtitel),
               numCats=2,catMethod=c(0:100),addLegend=FALSE,colourPalette=c("red","green"))

# ----------------

library(rworldmap)

options("scipen"=100, "digits"=4)

year <- 2015

target.year <- which(renewable.electricity.output$year==year)

head(renewable.electricity.output[target.year,])

renewable.electricity.output.year <- na.omit(renewable.electricity.output[target.year,])

renewable.electricity.output.year <- renewable.electricity.output.year[which(renewable.electricity.output.year$EG.ELC.RNEW.ZS>0),]

mapped_data <- joinCountryData2Map(renewable.electricity.output.year, joinCode = "ISO2", nameJoinColumn = "iso2c")

par(mai=c(0,0,0.2,0),xaxs="i",yaxs="i")

mapCountryData(mapped_data, nameColumnToPlot = "EG.ELC.RNEW.ZS",mapTitle=paste(titel,subtitel))

mapCountryData(mapped_data, nameColumnToPlot = "EG.ELC.RNEW.ZS",mapRegion="world",mapTitle=paste(titel,subtitel),

numCats=2,catMethod=c(0:100),addLegend=FALSE,colourPalette=c("red","green"))

mapCountryData(mapped_data, nameColumnToPlot = "EG.ELC.RNEW.ZS",mapRegion="world",mapTitle=paste(titel,subtitel),
               numCats=3,catMethod=c(0:100),addLegend=FALSE,colourPalette=c("red","yellow","green"))
legend(-179, -10, legend=c("low", "medium","high"), col=c("red", "yellow","green"), lty=1:3, cex=0.6)

mapCountryData(mapped_data, nameColumnToPlot = "EG.ELC.RNEW.ZS",mapRegion="world",mapTitle=paste(titel,subtitel),

numCats=3,catMethod=c(0:100),addLegend=FALSE,colourPalette=c("red","yellow","green"))

legend(-179, -10, legend=c("low", "medium","high"), col=c("red", "yellow","green"), lty=1:3, cex=0.6)

Sound System Optimisation Seminar – Certificate of Completion

02/01/2019AllgemeinMartin Stoppacher

access file systems – recursive function r

05/06/2017AllgemeinMartin Stoppacher

I was playing around with recursion. This is a short snippet about iterating through folder structures. It shows all directories and subdirectories within a file path.

library(rlang)  # is_empty() to check if a folder is empty

show.dir <- function(path="/",sleep.time=1){
  if(!is_empty(path)){
    
    print(path);
    print(dir(path,all.files = FALSE, full.names = FALSE, recursive = FALSE));
    Sys.sleep(sleep.time);
        
    #print(charToRaw(paste(dir(path,all.files = FALSE, full.names = FALSE, recursive = FALSE)))); # print HEX for fun

    current.dir <- dir(path,all.files = FALSE, full.names = TRUE, recursive = FALSE);
    
    for(i in 1:length(current.dir)){
      
      #print(i)
      #print(current.dir[i])
    
      if(!is_empty(dir(current.dir[i]))){
      
        #show.dir(current.dir[i])
        #print(dir(current.dir[i]));
        
        #print(charToRaw(paste(dir(current.dir[i]),collapse=""))) # print HEX for fun
        
        #Sys.sleep(sleep.time);
        show.dir(current.dir[i],sleep.time=sleep.time)
        
      }else{print(current.dir[i]);print("empty");}
      
    }
    
  }
}
show.dir(path="/Users/impac/Desktop/test",sleep.time=5)
show.dir()

show.dir(sleep.time=0.3)
show.dir(sleep.time=0.01)

library(rlang) # is_empty() to check if a folder is empty

show.dir <- function(path="/",sleep.time=1){

if(!is_empty(path)){

print(path);

print(dir(path,all.files = FALSE, full.names = FALSE, recursive = FALSE));

Sys.sleep(sleep.time);

#print(charToRaw(paste(dir(path,all.files = FALSE, full.names = FALSE, recursive = FALSE)))); # print HEX for fun

current.dir <- dir(path,all.files = FALSE, full.names = TRUE, recursive = FALSE);

for(i in 1:length(current.dir)){

#print(i)

#print(current.dir[i])

if(!is_empty(dir(current.dir[i]))){

#show.dir(current.dir[i])

#print(dir(current.dir[i]));

#print(charToRaw(paste(dir(current.dir[i]),collapse=""))) # print HEX for fun

#Sys.sleep(sleep.time);

show.dir(current.dir[i],sleep.time=sleep.time)

}else{print(current.dir[i]);print("empty");}

}

show.dir(path="/Users/impac/Desktop/test",sleep.time=5)

show.dir()

show.dir(sleep.time=0.3)

show.dir(sleep.time=0.01)

This is the output from a test directory:

[1] "/Users/impac/Desktop/test"
[1] "t1" "t2" "t3" "t4"
[1] "/Users/impac/Desktop/test/t1"
[1] "tt1" "tt2"
[1] "/Users/impac/Desktop/test/t1/tt1"
[1] "ttt1-3"
[1] "/Users/impac/Desktop/test/t1/tt1/ttt1-3"
[1] "empty"
[1] "/Users/impac/Desktop/test/t1/tt2"
[1] "ttt2-3"
[1] "/Users/impac/Desktop/test/t1/tt2/ttt2-3"
[1] "empty"
[1] "/Users/impac/Desktop/test/t2"
[1] "b1" "b2"
[1] "/Users/impac/Desktop/test/t2/b1"
[1] "empty"
[1] "/Users/impac/Desktop/test/t2/b2"
[1] "empty"
[1] "/Users/impac/Desktop/test/t3"
[1] "empty"
[1] "/Users/impac/Desktop/test/t4"
[1] "empty"

[1] "/Users/impac/Desktop/test"

[1] "t1" "t2" "t3" "t4"

[1] "/Users/impac/Desktop/test/t1"

[1] "tt1" "tt2"

[1] "/Users/impac/Desktop/test/t1/tt1"

[1] "ttt1-3"

[1] "/Users/impac/Desktop/test/t1/tt1/ttt1-3"

[1] "empty"

[1] "/Users/impac/Desktop/test/t1/tt2"

[1] "ttt2-3"

[1] "/Users/impac/Desktop/test/t1/tt2/ttt2-3"

[1] "empty"

[1] "/Users/impac/Desktop/test/t2"

[1] "b1" "b2"

[1] "/Users/impac/Desktop/test/t2/b1"

[1] "empty"

[1] "/Users/impac/Desktop/test/t2/b2"

[1] "empty"

[1] "/Users/impac/Desktop/test/t3"

[1] "empty"

[1] "/Users/impac/Desktop/test/t4"

[1] "empty"

Plane crash data

06/04/2016Allgemein, Data Analysis, R, StatisticsMartin Stoppacher

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# A Simple Skript for planecrash data from: 
# http://www.planecrashinfo.com/database.htm
# The aviation accident database includes: 
# All civil and commercial aviation accidents of scheduled and non-scheduled 
# passenger airliners worldwide, which resulted in a fatality 
# (including all U.S. Part 121 and Part 135 fatal accidents)
# All cargo, positioning, ferry and test flight fatal accidents.
# All military transport accidents with 10 or more fatalities.
# All commercial and military helicopter accidents with greater than 
# 10 fatalities.
# All civil and military airship accidents involving fatalities.
# Aviation accidents involving the death of famous people.
# Aviation accidents or incidents of noteworthy interest.

# Database Format
# Date:	 Date of accident,  in the format - January 01, 2001
# Time:	 Local time, in 24 hr. format unless otherwise specified
# Airline/Op:	 Airline or operator of the aircraft
# Flight #:	 Flight number assigned by the aircraft operator
# Route:	 Complete or partial route flown prior to the accident
# AC Type:	 Aircraft type
# Reg:	 ICAO registration of the aircraft
# cn / ln:	 Construction or serial number / Line or fuselage number
# Aboard:	 Total aboard (passengers / crew)
# Fatalities:	 Total fatalities aboard (passengers / crew)
# Ground:	 Total killed on the ground
# Summary:	 Brief description of the accident and cause if known


# to be done: histograms, t test?, types of chrashes, 

# country filter?, region filter, plane filter
# text analysis, 

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# packages

library("XML")
library("quantmod")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 

rm(list = ls(all = TRUE))
getwd()
#system("ls")
setwd("~/ownCloud/STA_Statistics/PlaneCrashData")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# get the available years

url <- "http://www.planecrashinfo.com/database.htm"
tables <- readHTMLTable(url,header=FALSE)
years <- NULL
for(i in 1:length(tables[[2]][,1])){
  for(e in 2:length(tables[[2]][1,])){
    years<-c(years,levels(tables[[2]][i,e])[i])
  }
}
years
years<-na.omit(as.numeric(years))
years

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# get all data tables

table.crash.all <- list()

for(i in 1:length(years)){
  url2 <- paste("http://www.planecrashinfo.com/",
                years[i],"/",years[i],".htm",sep="")
  assign(paste("table.crash.",years[i],sep=""),readHTMLTable(url2))
  table.crash.all[[i]] <- get(paste("table.crash.",years[i],sep=""))
  print(paste("Year",years[i],"downloaded"))
  }
#table.crash.all
length(table.crash.all)
save(table.crash.all,file="planecrashdata_list.RData")
for(i in 1:length(years)){
  data<-table.crash.all[[i]]
  save(data,file=paste("planecrashdata_",years[i],".RData",sep=""))
}
for(i in 1:length(years)){
  data<-table.crash.all[[i]]
  write.csv(data,file=paste("planecrashdata_",years[i],".csv",sep=""))
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Plane Crashes per Year

years.crashes <- data.frame()
for(i in 1:length(years)){
  years.crashes[i,1] <- years[i]
  years.crashes[i,2] <- length(as.data.frame(table.crash.all[[i]])[,1])
}
years.crashes
rownames(years.crashes)<-years.crashes[,1]

par(mfrow=c(1,1))

plot(years.crashes,type="b",main="Plane Crashes per Year",
     xlab="Year",ylab="number of crashes")
points(years.crashes,col="red")

par(mfrow=c(1,1))

setwd("/Users/impac/ownCloud/STA_Statistics/PlaneCrashData/plots")

jpeg(file="Plane Crashes per Year-points.jpeg", width = 1080, 
     height = 1080, pointsize = 12, quality = 75)

plot(years.crashes,type="b",main="Plane Crashes per Year",
     xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")
dev.off()

jpeg(file="Plane Crashes per Year-bars.jpeg", width = 1080, 
     height = 1080, pointsize = 12, quality = 75)

plot(years.crashes,type="h",main="Plane Crashes per Year",
     xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")
dev.off()

jpeg(file="Plane Crashes per Year-steps.jpeg", width = 1080, 
     height = 1080, pointsize = 12, quality = 75)

plot(years.crashes,type="s",main="Plane Crashes per Year",
     xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")
dev.off()

par(mfrow=c(2,2))

plot(years.crashes[1:30,],type="h",main=paste("Plane Crashes per Year",
                           years.crashes[1,1],"-",years.crashes[30,1]),xlab="Year",
                           ylab="number of crashes")

points(years.crashes[1:30,],col="red")
abline(h=mean(years.crashes[1:30,2]),col="blue")

plot(years.crashes[31:60,],type="h",main=paste("Plane Crashes per Year",years.crashes[31,1],"-",years.crashes[60,1]),xlab="Year",ylab="number of crashes")
points(years.crashes[31:60,],col="red")
abline(h=mean(years.crashes[31:60,2]),col="blue")

plot(years.crashes[61:80,],type="h",main=paste("Plane Crashes per Year",years.crashes[61,1],"-",years.crashes[80,1]),xlab="Year",ylab="number of crashes")
points(years.crashes[61:80,],col="red")
abline(h=mean(years.crashes[61:80,2]),col="blue")

plot(years.crashes[81:length(years.crashes[,1]),],type="h",main=paste("Plane Crashes per Year",years.crashes[81,1],"-",years.crashes[length(years.crashes[,1]),1]),xlab="Year",ylab="number of crashes")
points(years.crashes[81:length(years.crashes[,1]),],col="red")
abline(h=mean(years.crashes[81:length(years.crashes[,1]),2]),col="blue")


#jpeg(file="Plane Crashes per Year.jpeg", width = 1080, height = 1080, pointsize = 12, quality = 75)
#par(mfrow=c(2,2))

#plot(years.crashes[1:30,],type="h",main=paste("Plane Crashes per Year",years.crashes[1,1],"-",years.crashes[30,1]),xlab="Year",ylab="number of crashes")
#points(years.crashes[1:30,],col="red")
#abline(h=mean(years.crashes[1:30,2]),col="blue")

#plot(years.crashes[31:60,],type="h",main=paste("Plane Crashes per Year",years.crashes[31,1],"-",years.crashes[60,1]),xlab="Year",ylab="number of crashes")
#points(years.crashes[31:60,],col="red")
#abline(h=mean(years.crashes[31:60,2]),col="blue")

#plot(years.crashes[61:80,],type="h",main=paste("Plane Crashes per Year",years.crashes[61,1],"-",years.crashes[80,1]),xlab="Year",ylab="number of crashes")
#points(years.crashes[61:80,],col="red")
#abline(h=mean(years.crashes[61:80,2]),col="blue")

#plot(years.crashes[81:length(years.crashes[,1]),],type="h",main=paste("Plane Crashes per Year",years.crashes[81,1],"-",years.crashes[length(years.crashes[,1]),1]),xlab="Year",ylab="number of crashes")
#points(years.crashes[81:length(years.crashes[,1]),],col="red")
#abline(h=mean(years.crashes[81:length(years.crashes[,1]),2]),col="blue")

#dev.off()

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Histograms of the number of crashes

par(mfrow=c(1,1))

hist(years.crashes[,2], nclass=length(years.crashes[,2]), 
     freq=TRUE, main='histogram of plane crashes',
     xlab="number of chrashes per year", col = "lightblue", border = "blue")

hist(years.crashes[,2], nclass=length(years.crashes[,2]), freq=FALSE, 
     main='histogram of plane crashes',xlab="number of chrashes per year", 
     col = "lightblue", border = "blue")

curve(dnorm(x, mean=mean(years.crashes[,2]),sd=sd(years.crashes[,2])), 
      add=TRUE, col="red")

par(mfrow=c(1,2))

qqnorm(years.crashes[,2])
qqline(years.crashes[,2], col = 2)
plot(density(years.crashes[,2]),log='y', main='')
curve(dnorm(x, mean=mean(years.crashes[,2]),
            sd=sd(years.crashes[,2])), log="y", add=TRUE, col="red")

# different periods

par(mfrow=c(1,2))

# 1980-1999
#hist(years.crashes[61:80,2], nclass=20, freq=TRUE, 
#     main='',xlab="number of chrashes per year", col = "lightblue", 
#     border = "blue")

hist(years.crashes[61:80,2], nclass=20, freq=FALSE, main='1980-1999',xlab="number of chrashes per year", col = "lightblue", border = "blue")
curve(dnorm(x, mean=mean(years.crashes[61:80,2]),sd=sd(years.crashes[61:80,2])), add=TRUE, col="red")

# 2000-2020
#hist(years.crashes[81:length(years.crashes[,1]),2], nclass=length(years.crashes[81:length(years.crashes[,1]),2]), freq=TRUE, main='',xlab="number of chrashes per year", col = "lightblue", border = "blue")
hist(years.crashes[81:length(years.crashes[,1]),2], nclass=length(years.crashes[81:length(years.crashes[,1]),2]), freq=FALSE, main='2000-2020',xlab="number of chrashes per year", col = "lightblue", border = "blue")
curve(dnorm(x,mean=mean(years.crashes[81:length(years.crashes[,1]),2]),sd=sd(years.crashes[81:length(years.crashes[,1]),2])), add=TRUE, col="red")

par(mfrow=c(2,2))

qqnorm(years.crashes[61:80,2])
qqline(years.crashes[61:80,2], col = 2)
plot(density(years.crashes[61:80,2]),log='y', main='')
curve(dnorm(x, mean=mean(years.crashes[61:80,2]),sd=sd(years.crashes[61:80,2])), log="y", add=TRUE, col="red")

qqnorm(years.crashes[81:length(years.crashes[,1]),2])
qqline(years.crashes[81:length(years.crashes[,1]),2], col = 2)
plot(density(years.crashes[81:length(years.crashes[,1]),2]),log='y', main='')
curve(dnorm(x, mean=mean(years.crashes[81:length(years.crashes[,1]),2]),sd=sd(years.crashes[81:length(years.crashes[,1]),2])), log="y", add=TRUE, col="red")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# boxplots

par(mfrow=c(1,1))

boxplot(years.crashes[,2], main="number of plain crashes per year 1920 - 2020", xlab="", ylab="plain crashes per year")

boxplot(years.crashes[,2], notch=TRUE, col=(c("gold","darkgreen")), main="number of plain crashes per year 1920 - 2020", ylab="plain crashes per year")

years.crashes[1:30,3]<-rep("1920-1949",30)
years.crashes[31:60,3]<-rep("1950-1979",30)278750
years.crashes[61:80,3]<-rep("1980-1999",20)
years.crashes[81:length(years.crashes[,1]),3]<-rep("2000-2020",21)

colnames(years.crashes)<-c("year","ncrashes","cat")

boxplot(ncrashes~cat, data=years.crashes, notch=TRUE, col=(c("green","blue","blue","green")), main="number of plain crashes per year - groups", xlab="")

# start from here
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# t test

# independent 2-group t-test    t.test(y~x) # where y is numeric and x is a binary factor
# independent 2-group t-test    t.test(y1,y2) # where y1 and y2 are numeric
# paired t-test                 t.test(y1,y2,paired=TRUE) # where y1 & y2 are numeric
# one sample t-test             t.test(y,mu=3) # Ho: mu=3

t.test(years.crashes[61:80,2],years.crashes[81:length(years.crashes[,1]),2])

# Welch Two Sample t-test
#
# data:  years.crashes[61:80, 2] and years.crashes[81:length(years.crashes[, 1]), 2]
# t = 5.4142, df = 31.453, p-value = 6.268e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   14.70178 32.45536
# sample estimates:
#   mean of x mean of y 
# 59.15000  35.57143 

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# rate of difference

par(mfrow=c(3,1))

rate<-na.omit(ROC(years.crashes[,2]))
plot(rate,type="l")
plot(cumsum(rate),type="l")
hist(rate, nclass=20)

rate<-na.omit(ROC(years.crashes[61:80,2]))
plot(rate,type="l")
plot(cumsum(rate),type="l")
hist(rate, nclass=20)

rate<-na.omit(ROC(years.crashes[81:length(years.crashes[,1]),2]))
plot(rate,type="l")
plot(cumsum(rate),type="l")
hist(rate, nclass=20)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# One big Table

table.crash.one.table <- data.frame(table.crash.all[[1]],stringsAsFactors = TRUE)
colnames(table.crash.one.table)<-c("date","location","type","fatalities")

for(i in 2:length(years)){
  data <- data.frame(table.crash.all[[i]],stringsAsFactors = TRUE)
  colnames(data)<-c("date","location","type","fatalities")
  table.crash.one.table <- rbind(table.crash.one.table,data)
}
table.crash.one.table

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# create a time series object

date <- gsub("Jan","01",table.crash.one.table[,1])
date <- gsub("Feb","02",date)
date <- gsub("Mar","03",date)
date <- gsub("Apr","04",date)
date <- gsub("May","05",date)
date <- gsub("Jun","06",date)
date <- gsub("Jul","07",date)
date <- gsub("Aug","08",date)
date <- gsub("Sep","09",date)
date <- gsub("Oct","10",date)
date <- gsub("Nov","11",date)
date <- gsub("Dec","12",date)

gsub("/","-",table.crash.one.table[,4])

fatalities <- strsplit(gsub("/","-",table.crash.one.table[,4]),'-')
number<-NULL
for(i in 1:length(fatalities)){number<-c(number,fatalities[[i]][1])}
#as.numeric(number)
table.crash.one.table[,5]<-as.numeric(number)


#strsplit(gsub("(","-",fatalities[[i]][2], fixed = TRUE),"-")[1]
number2<-NULL
for(i in 1:length(fatalities)){number2<-c(number2,
                strsplit(gsub("(","-",fatalities[[i]][2], fixed = TRUE),"-")[[1]][1])}
number2
table.crash.one.table[,6]<-as.numeric(number2)

colnames(table.crash.one.table)<-c("date","location","type",
                                   "fatalities","fatalities-passengers","fatalities-crew")

#table.crash.one.table.xts <- xts(table.crash.one.table[,2:6], order.by=as.Date(date,"%d %m %Y"))
#colnames(table.crash.one.table.xts)<-c("location","type","fatalities",
#"passengerfatalities","crewfatalities")

table.crash.one.table.xts <- xts(table.crash.one.table[,5:6], order.by=as.Date(date,"%d %m %Y"))
colnames(table.crash.one.table.xts)<-c("totalfatalities","totalaboard")

table.crash.one.table.xts$totalfatalities
table.crash.one.table.xts$totalaboard

plot(table.crash.one.table.xts$totalfatalities, main="Plane Crash Total Fatalities")

plot(table.crash.one.table.xts$totalaboard, main="Plane Crash Total Aboard")

plot(table.crash.one.table.xts$totalaboard, main="Plane Crash Total Aboard vs Total Fatalities")
lines(table.crash.one.table.xts$totalfatalities,col="red")

plot(table.crash.one.table.xts$totalaboard["2000/2020"], 
     main="Plane Crash Total Aboard vs Total Fatalities",type="c")

plot(table.crash.one.table.xts$totalaboard["2000/2020"], 
     main="Plane Crash Total Aboard vs Total Fatalities",type="o")

plot(table.crash.one.table.xts$totalaboard["2000/2020"], 
     main="Plane Crash Total Aboard vs Total Fatalities",type="h")

plot(table.crash.one.table.xts$totalaboard["2000/2020"], 
     main="Plane Crash Total Aboard vs Total Fatalities",type="b")
points(table.crash.one.table.xts$totalfatalities["2000/2020"],col="red")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Total killed on the ground



number3<-NULL
for(i in 1:length(fatalities)){
  sub <- gsub("(","-",fatalities[[i]][2], fixed = TRUE)
  sub <- gsub(")","-",sub, fixed = TRUE)
  number3<-c(number3,strsplit(sub,"-")[[1]][2])
}
number3
table.crash.one.table[,7]<-as.numeric(number3)

colnames(table.crash.one.table)<-c("date","location","type",
                          "fatalities","fatalities-passengers","fatalities-crew","totalkilledontheground")

table.crash.one.table.xts <- xts(table.crash.one.table[,5:7], order.by=as.Date(date,"%d %m %Y"))
colnames(table.crash.one.table.xts)<-c("totalfatalities","totalaboard","totalkilledontheground")


plot(table.crash.one.table.xts$totalkilledontheground, main="Plane Crash Total killed on the ground")
plot(table.crash.one.table.xts$totalkilledontheground, main="Plane Crash Total killed on the ground",type="b")
plot(table.crash.one.table.xts$totalkilledontheground, main="Plane Crash Total killed on the ground",type="c")

# Martin Stoppacher                                                             #
# office@martinstoppacher.com                                                   #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#################################################################################

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# A Simple Skript for planecrash data from:

# http://www.planecrashinfo.com/database.htm

# The aviation accident database includes:

# All civil and commercial aviation accidents of scheduled and non-scheduled

# passenger airliners worldwide, which resulted in a fatality

# (including all U.S. Part 121 and Part 135 fatal accidents)

# All cargo, positioning, ferry and test flight fatal accidents.

# All military transport accidents with 10 or more fatalities.

# All commercial and military helicopter accidents with greater than

# 10 fatalities.

# All civil and military airship accidents involving fatalities.

# Aviation accidents involving the death of famous people.

# Aviation accidents or incidents of noteworthy interest.

# Database Format

# Date: Date of accident, in the format - January 01, 2001

# Time: Local time, in 24 hr. format unless otherwise specified

# Airline/Op: Airline or operator of the aircraft

# Flight #: Flight number assigned by the aircraft operator

# Route: Complete or partial route flown prior to the accident

# AC Type: Aircraft type

# Reg: ICAO registration of the aircraft

# cn / ln: Construction or serial number / Line or fuselage number

# Aboard: Total aboard (passengers / crew)

# Fatalities: Total fatalities aboard (passengers / crew)

# Ground: Total killed on the ground

# Summary: Brief description of the accident and cause if known

# to be done: histograms, t test?, types of chrashes,

# country filter?, region filter, plane filter

# text analysis,

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# packages

library("XML")

library("quantmod")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

rm(list = ls(all = TRUE))

getwd()

#system("ls")

setwd("~/ownCloud/STA_Statistics/PlaneCrashData")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# get the available years

url <- "http://www.planecrashinfo.com/database.htm"

tables <- readHTMLTable(url,header=FALSE)

years <- NULL

for(i in 1:length(tables[[2]][,1])){

for(e in 2:length(tables[[2]][1,])){

years<-c(years,levels(tables[[2]][i,e])[i])

}

years

years<-na.omit(as.numeric(years))

years

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# get all data tables

table.crash.all <- list()

for(i in 1:length(years)){

url2 <- paste("http://www.planecrashinfo.com/",

years[i],"/",years[i],".htm",sep="")

assign(paste("table.crash.",years[i],sep=""),readHTMLTable(url2))

table.crash.all[[i]] <- get(paste("table.crash.",years[i],sep=""))

print(paste("Year",years[i],"downloaded"))

}

#table.crash.all

length(table.crash.all)

save(table.crash.all,file="planecrashdata_list.RData")

for(i in 1:length(years)){

data<-table.crash.all[[i]]

save(data,file=paste("planecrashdata_",years[i],".RData",sep=""))

}

for(i in 1:length(years)){

data<-table.crash.all[[i]]

write.csv(data,file=paste("planecrashdata_",years[i],".csv",sep=""))

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Plane Crashes per Year

years.crashes <- data.frame()

for(i in 1:length(years)){

years.crashes[i,1] <- years[i]

years.crashes[i,2] <- length(as.data.frame(table.crash.all[[i]])[,1])

}

years.crashes

rownames(years.crashes)<-years.crashes[,1]

par(mfrow=c(1,1))

plot(years.crashes,type="b",main="Plane Crashes per Year",

xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")

par(mfrow=c(1,1))

setwd("/Users/impac/ownCloud/STA_Statistics/PlaneCrashData/plots")

jpeg(file="Plane Crashes per Year-points.jpeg", width = 1080,

height = 1080, pointsize = 12, quality = 75)

plot(years.crashes,type="b",main="Plane Crashes per Year",

xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")

dev.off()

jpeg(file="Plane Crashes per Year-bars.jpeg", width = 1080,

height = 1080, pointsize = 12, quality = 75)

plot(years.crashes,type="h",main="Plane Crashes per Year",

xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")

dev.off()

jpeg(file="Plane Crashes per Year-steps.jpeg", width = 1080,

height = 1080, pointsize = 12, quality = 75)

plot(years.crashes,type="s",main="Plane Crashes per Year",

xlab="Year",ylab="number of crashes")

points(years.crashes,col="red")

dev.off()

par(mfrow=c(2,2))

plot(years.crashes[1:30,],type="h",main=paste("Plane Crashes per Year",

years.crashes[1,1],"-",years.crashes[30,1]),xlab="Year",

ylab="number of crashes")

points(years.crashes[1:30,],col="red")

abline(h=mean(years.crashes[1:30,2]),col="blue")

plot(years.crashes[31:60,],type="h",main=paste("Plane Crashes per Year",years.crashes[31,1],"-",years.crashes[60,1]),xlab="Year",ylab="number of crashes")

points(years.crashes[31:60,],col="red")

abline(h=mean(years.crashes[31:60,2]),col="blue")

plot(years.crashes[61:80,],type="h",main=paste("Plane Crashes per Year",years.crashes[61,1],"-",years.crashes[80,1]),xlab="Year",ylab="number of crashes")

points(years.crashes[61:80,],col="red")

abline(h=mean(years.crashes[61:80,2]),col="blue")

plot(years.crashes[81:length(years.crashes[,1]),],type="h",main=paste("Plane Crashes per Year",years.crashes[81,1],"-",years.crashes[length(years.crashes[,1]),1]),xlab="Year",ylab="number of crashes")

points(years.crashes[81:length(years.crashes[,1]),],col="red")

abline(h=mean(years.crashes[81:length(years.crashes[,1]),2]),col="blue")

#jpeg(file="Plane Crashes per Year.jpeg", width = 1080, height = 1080, pointsize = 12, quality = 75)

#par(mfrow=c(2,2))

#plot(years.crashes[1:30,],type="h",main=paste("Plane Crashes per Year",years.crashes[1,1],"-",years.crashes[30,1]),xlab="Year",ylab="number of crashes")

#points(years.crashes[1:30,],col="red")

#abline(h=mean(years.crashes[1:30,2]),col="blue")

#plot(years.crashes[31:60,],type="h",main=paste("Plane Crashes per Year",years.crashes[31,1],"-",years.crashes[60,1]),xlab="Year",ylab="number of crashes")

#points(years.crashes[31:60,],col="red")

#abline(h=mean(years.crashes[31:60,2]),col="blue")

#plot(years.crashes[61:80,],type="h",main=paste("Plane Crashes per Year",years.crashes[61,1],"-",years.crashes[80,1]),xlab="Year",ylab="number of crashes")

#points(years.crashes[61:80,],col="red")

#abline(h=mean(years.crashes[61:80,2]),col="blue")

#plot(years.crashes[81:length(years.crashes[,1]),],type="h",main=paste("Plane Crashes per Year",years.crashes[81,1],"-",years.crashes[length(years.crashes[,1]),1]),xlab="Year",ylab="number of crashes")

#points(years.crashes[81:length(years.crashes[,1]),],col="red")

#abline(h=mean(years.crashes[81:length(years.crashes[,1]),2]),col="blue")

#dev.off()

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Histograms of the number of crashes

par(mfrow=c(1,1))

hist(years.crashes[,2], nclass=length(years.crashes[,2]),

freq=TRUE, main='histogram of plane crashes',

xlab="number of chrashes per year", col = "lightblue", border = "blue")

hist(years.crashes[,2], nclass=length(years.crashes[,2]), freq=FALSE,

main='histogram of plane crashes',xlab="number of chrashes per year",

col = "lightblue", border = "blue")

curve(dnorm(x, mean=mean(years.crashes[,2]),sd=sd(years.crashes[,2])),

add=TRUE, col="red")

par(mfrow=c(1,2))

qqnorm(years.crashes[,2])

qqline(years.crashes[,2], col = 2)

plot(density(years.crashes[,2]),log='y', main='')

curve(dnorm(x, mean=mean(years.crashes[,2]),

sd=sd(years.crashes[,2])), log="y", add=TRUE, col="red")

# different periods

par(mfrow=c(1,2))

# 1980-1999

#hist(years.crashes[61:80,2], nclass=20, freq=TRUE,

# main='',xlab="number of chrashes per year", col = "lightblue",

# border = "blue")

hist(years.crashes[61:80,2], nclass=20, freq=FALSE, main='1980-1999',xlab="number of chrashes per year", col = "lightblue", border = "blue")

curve(dnorm(x, mean=mean(years.crashes[61:80,2]),sd=sd(years.crashes[61:80,2])), add=TRUE, col="red")

# 2000-2020

#hist(years.crashes[81:length(years.crashes[,1]),2], nclass=length(years.crashes[81:length(years.crashes[,1]),2]), freq=TRUE, main='',xlab="number of chrashes per year", col = "lightblue", border = "blue")

hist(years.crashes[81:length(years.crashes[,1]),2], nclass=length(years.crashes[81:length(years.crashes[,1]),2]), freq=FALSE, main='2000-2020',xlab="number of chrashes per year", col = "lightblue", border = "blue")

curve(dnorm(x,mean=mean(years.crashes[81:length(years.crashes[,1]),2]),sd=sd(years.crashes[81:length(years.crashes[,1]),2])), add=TRUE, col="red")

par(mfrow=c(2,2))

qqnorm(years.crashes[61:80,2])

qqline(years.crashes[61:80,2], col = 2)

plot(density(years.crashes[61:80,2]),log='y', main='')

curve(dnorm(x, mean=mean(years.crashes[61:80,2]),sd=sd(years.crashes[61:80,2])), log="y", add=TRUE, col="red")

qqnorm(years.crashes[81:length(years.crashes[,1]),2])

qqline(years.crashes[81:length(years.crashes[,1]),2], col = 2)

plot(density(years.crashes[81:length(years.crashes[,1]),2]),log='y', main='')

curve(dnorm(x, mean=mean(years.crashes[81:length(years.crashes[,1]),2]),sd=sd(years.crashes[81:length(years.crashes[,1]),2])), log="y", add=TRUE, col="red")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# boxplots

par(mfrow=c(1,1))

boxplot(years.crashes[,2], main="number of plain crashes per year 1920 - 2020", xlab="", ylab="plain crashes per year")

boxplot(years.crashes[,2], notch=TRUE, col=(c("gold","darkgreen")), main="number of plain crashes per year 1920 - 2020", ylab="plain crashes per year")

years.crashes[1:30,3]<-rep("1920-1949",30)

years.crashes[31:60,3]<-rep("1950-1979",30)278750

years.crashes[61:80,3]<-rep("1980-1999",20)

years.crashes[81:length(years.crashes[,1]),3]<-rep("2000-2020",21)

colnames(years.crashes)<-c("year","ncrashes","cat")

boxplot(ncrashes~cat, data=years.crashes, notch=TRUE, col=(c("green","blue","blue","green")), main="number of plain crashes per year - groups", xlab="")

# start from here

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# t test

# independent 2-group t-test t.test(y~x) # where y is numeric and x is a binary factor

# independent 2-group t-test t.test(y1,y2) # where y1 and y2 are numeric

# paired t-test t.test(y1,y2,paired=TRUE) # where y1 & y2 are numeric

# one sample t-test t.test(y,mu=3) # Ho: mu=3

t.test(years.crashes[61:80,2],years.crashes[81:length(years.crashes[,1]),2])

# Welch Two Sample t-test

# data: years.crashes[61:80, 2] and years.crashes[81:length(years.crashes[, 1]), 2]

# t = 5.4142, df = 31.453, p-value = 6.268e-06

# alternative hypothesis: true difference in means is not equal to 0

# 95 percent confidence interval:

# 14.70178 32.45536

# sample estimates:

# mean of x mean of y

# 59.15000 35.57143

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# rate of difference

par(mfrow=c(3,1))

rate<-na.omit(ROC(years.crashes[,2]))

plot(rate,type="l")

plot(cumsum(rate),type="l")

hist(rate, nclass=20)

rate<-na.omit(ROC(years.crashes[61:80,2]))

plot(rate,type="l")

plot(cumsum(rate),type="l")

hist(rate, nclass=20)

rate<-na.omit(ROC(years.crashes[81:length(years.crashes[,1]),2]))

plot(rate,type="l")

plot(cumsum(rate),type="l")

hist(rate, nclass=20)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# One big Table

table.crash.one.table <- data.frame(table.crash.all[[1]],stringsAsFactors = TRUE)

colnames(table.crash.one.table)<-c("date","location","type","fatalities")

for(i in 2:length(years)){

data <- data.frame(table.crash.all[[i]],stringsAsFactors = TRUE)

colnames(data)<-c("date","location","type","fatalities")

table.crash.one.table <- rbind(table.crash.one.table,data)

}

table.crash.one.table

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# create a time series object

date <- gsub("Jan","01",table.crash.one.table[,1])

date <- gsub("Feb","02",date)

date <- gsub("Mar","03",date)

date <- gsub("Apr","04",date)

date <- gsub("May","05",date)

date <- gsub("Jun","06",date)

date <- gsub("Jul","07",date)

date <- gsub("Aug","08",date)

date <- gsub("Sep","09",date)

date <- gsub("Oct","10",date)

date <- gsub("Nov","11",date)

date <- gsub("Dec","12",date)

gsub("/","-",table.crash.one.table[,4])

fatalities <- strsplit(gsub("/","-",table.crash.one.table[,4]),'-')

number<-NULL

for(i in 1:length(fatalities)){number<-c(number,fatalities[[i]][1])}

#as.numeric(number)

table.crash.one.table[,5]<-as.numeric(number)

#strsplit(gsub("(","-",fatalities[[i]][2], fixed = TRUE),"-")[1]

number2<-NULL

for(i in 1:length(fatalities)){number2<-c(number2,

strsplit(gsub("(","-",fatalities[[i]][2], fixed = TRUE),"-")[[1]][1])}

number2

table.crash.one.table[,6]<-as.numeric(number2)

colnames(table.crash.one.table)<-c("date","location","type",

"fatalities","fatalities-passengers","fatalities-crew")

#table.crash.one.table.xts <- xts(table.crash.one.table[,2:6], order.by=as.Date(date,"%d %m %Y"))

#colnames(table.crash.one.table.xts)<-c("location","type","fatalities",

#"passengerfatalities","crewfatalities")

table.crash.one.table.xts <- xts(table.crash.one.table[,5:6], order.by=as.Date(date,"%d %m %Y"))

colnames(table.crash.one.table.xts)<-c("totalfatalities","totalaboard")

table.crash.one.table.xts$totalfatalities

table.crash.one.table.xts$totalaboard

plot(table.crash.one.table.xts$totalfatalities, main="Plane Crash Total Fatalities")

plot(table.crash.one.table.xts$totalaboard, main="Plane Crash Total Aboard")

plot(table.crash.one.table.xts$totalaboard, main="Plane Crash Total Aboard vs Total Fatalities")

lines(table.crash.one.table.xts$totalfatalities,col="red")

plot(table.crash.one.table.xts$totalaboard["2000/2020"],

main="Plane Crash Total Aboard vs Total Fatalities",type="c")

plot(table.crash.one.table.xts$totalaboard["2000/2020"],

main="Plane Crash Total Aboard vs Total Fatalities",type="o")

plot(table.crash.one.table.xts$totalaboard["2000/2020"],

main="Plane Crash Total Aboard vs Total Fatalities",type="h")

plot(table.crash.one.table.xts$totalaboard["2000/2020"],

main="Plane Crash Total Aboard vs Total Fatalities",type="b")

points(table.crash.one.table.xts$totalfatalities["2000/2020"],col="red")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Total killed on the ground

number3<-NULL

for(i in 1:length(fatalities)){

sub <- gsub("(","-",fatalities[[i]][2], fixed = TRUE)

sub <- gsub(")","-",sub, fixed = TRUE)

number3<-c(number3,strsplit(sub,"-")[[1]][2])

}

number3

table.crash.one.table[,7]<-as.numeric(number3)

colnames(table.crash.one.table)<-c("date","location","type",

"fatalities","fatalities-passengers","fatalities-crew","totalkilledontheground")

table.crash.one.table.xts <- xts(table.crash.one.table[,5:7], order.by=as.Date(date,"%d %m %Y"))

colnames(table.crash.one.table.xts)<-c("totalfatalities","totalaboard","totalkilledontheground")

plot(table.crash.one.table.xts$totalkilledontheground, main="Plane Crash Total killed on the ground")

plot(table.crash.one.table.xts$totalkilledontheground, main="Plane Crash Total killed on the ground",type="b")

plot(table.crash.one.table.xts$totalkilledontheground, main="Plane Crash Total killed on the ground",type="c")

# Martin Stoppacher #

# office@martinstoppacher.com #

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

#################################################################################

My Sounds (Freesound.org)

01/31/2016Allgemein, Engineering, SoundengineeringMartin Stoppacher

This is a collection of my sounds at freesound.org. I try to keep the list up to date but sometimes there are more sounds available on freesound therefore please check out my freesound account and freesound in general. “Freesound is a collaborative database of Creative Commons Licensed sounds. Browse, download and share sounds.” – freesound.org

Support this Blog

bifurcation tree with R

06/23/2015Allgemein, Programming, R, StatisticsMartin Stoppacher

Lately i was exploring the logistic map function because i was fascinated again how chaotic behavior can arise from very simple circumstances (i.e. a rather simple equation).
Of course, i learned about this in school but never got a chance to write some code.
So here it is. 🙂

And here is some of the code i wrote:

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# bifurcation tree, Logistic map
# xn1=rxn(1-xn)

rm(list = ls(all = TRUE))

r <- 2.2           # rate
x <- 0.5           # population

(r*x*(1-x))

r <- 2.2
(x <- seq(100,1)/100)
y = r*x*(1-x)
plot(y)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# bifurcation tree, Logistic map

# xn1=rxn(1-xn)

rm(list = ls(all = TRUE))

r <- 2.2 # rate

x <- 0.5 # population

(r*x*(1-x))

r <- 2.2

(x <- seq(100,1)/100)

y = r*x*(1-x)

plot(y)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# if the starting population is 0.1

r <- 2.6
x <- 0.1
for(i in 2:20){
  (x[i] <- r*x[i-1]*(1-x[i-1]))
  print(i)
  print(x[i])
  #Sys.sleep(0.5)
}
x
plot(x, type="b")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# if the starting population is 0.1

r <- 2.6

x <- 0.1

for(i in 2:20){

(x[i] <- r*x[i-1]*(1-x[i-1]))

print(i)

print(x[i])

#Sys.sleep(0.5)

}

plot(x, type="b")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# if the starting population is 0.8

r <- 2.6
x <- 0.8
for(i in 2:20){
  (x[i] <- r*x[i-1]*(1-x[i-1]))
  print(i)
  print(x[i])
  #Sys.sleep(0.5)
}
x
plot(x, type="b")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# if the starting population is 0.8

r <- 2.6

x <- 0.8

for(i in 2:20){

(x[i] <- r*x[i-1]*(1-x[i-1]))

print(i)

print(x[i])

#Sys.sleep(0.5)

}

plot(x, type="b")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# if we make variations to the initial population all of them stabilize

plot(1:20,(1:20)/20,ylab="",xlab="",cex=0)

for(e in ((1:9)/10)){
r <- 2.6
x <- e
for(i in 2:20){
  (x[i] <- r*x[i-1]*(1-x[i-1]))
  print(i)
  print(x[i])
  #Sys.sleep(0.5)
}
x
lines(x, type="b")
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# if we make variations to the initial population all of them stabilize

plot(1:20,(1:20)/20,ylab="",xlab="",cex=0)

for(e in ((1:9)/10)){

r <- 2.6

x <- e

for(i in 2:20){

(x[i] <- r*x[i-1]*(1-x[i-1]))

print(i)

print(x[i])

#Sys.sleep(0.5)

}

lines(x, type="b")

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# what does this look like if there are variations to the groth rate?

plot(0,0, ylab="", xlab="", cex=0, xlim=c(1,20), ylim=c(0,1))

for(e in (1:30)/10){
  r <- e
  x <- 0.5
  for(i in 2:20){
    (x[i] <- r*x[i-1]*(1-x[i-1]))
    print(i)
    print(x[i])
    #Sys.sleep(0.5)
  }
  x
  lines(x, type="b")
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# what does this look like if there are variations to the groth rate?

plot(0,0, ylab="", xlab="", cex=0, xlim=c(1,20), ylim=c(0,1))

for(e in (1:30)/10){

r <- e

x <- 0.5

for(i in 2:20){

(x[i] <- r*x[i-1]*(1-x[i-1]))

print(i)

print(x[i])

#Sys.sleep(0.5)

}

lines(x, type="b")

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# ploting more values but only the last

plot(0,0, ylab="", xlab="", cex=0, xlim=c(2.5,4), ylim=c(0,1))

n<-300
m<-30
for(e in ((290:400)/100)){
  print(r)
  r <- e
  x <- 0.1
  for(i in 2:n){
    (x[i] <- r*x[i-1]*(1-x[i-1]))
    print(r)
    print(x[i])
  }
  x
  points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2)
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# ploting more values but only the last

plot(0,0, ylab="", xlab="", cex=0, xlim=c(2.5,4), ylim=c(0,1))

n<-300

m<-30

for(e in ((290:400)/100)){

print(r)

r <- e

x <- 0.1

for(i in 2:n){

(x[i] <- r*x[i-1]*(1-x[i-1]))

print(r)

print(x[i])

}

points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2)

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# adding color

plot(0,0, ylab="", xlab="", cex=0, xlim=c(0,4), ylim=c(0,1))

n<-100
m<-30
for(e in ((1:400)/100)){
  print(r)
  r <- e
  x <- 0.6
  for(i in 2:n){
    (x[i] <- r*x[i-1]*(1-x[i-1]))
    print(r)
    print(x[i])
  }
  x
  if(r<1)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="red")
  if(r>=1&r<2)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="black")
  if(r>=2&r<5)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="orange")
  if(r>=3&r<=3.5)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="blue")
  if(r>=3.5&r<=4)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="green")
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# adding color

plot(0,0, ylab="", xlab="", cex=0, xlim=c(0,4), ylim=c(0,1))

n<-100

m<-30

for(e in ((1:400)/100)){

print(r)

r <- e

x <- 0.6

for(i in 2:n){

(x[i] <- r*x[i-1]*(1-x[i-1]))

print(r)

print(x[i])

}

if(r<1)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="red")

if(r>=1&r<2)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="black")

if(r>=2&r<5)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="orange")

if(r>=3&r<=3.5)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="blue")

if(r>=3.5&r<=4)points(rep(r,(m+1)),x[(n-m):n], type="p", cex=0.2, lwd=0.2, col="green")

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# bifurcation tree as a function

log.tree <- function(r=2, x=0.4, n=200){
  for(i in 2:(n+1)){
    (x[i] <- r*x[i-1]*(1-x[i-1]))
  }
  x <- x[2:(n+1)]
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# bifurcation tree as a function

log.tree <- function(r=2, x=0.4, n=200){

for(i in 2:(n+1)){

(x[i] <- r*x[i-1]*(1-x[i-1]))

}

x <- x[2:(n+1)]

}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# using the function with sapply

point.count <- 50    # sets how many of the last observations to use
r.seq <- seq(1,4,by=0.001) # a sequence for different rates
plot(0,0, ylab="", xlab="", cex=0, xlim=c(0,4), ylim=c(0,1))
log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=3000),point.count)
for(i in 1:length(r.seq)) points(rep(r.seq[i],point.count),log.tree.points[,i], cex=0.2, lwd=0.2)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 

system.time(log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=300),point.count))
system.time(log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=3000),point.count))
system.time(log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=30000),point.count))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# using the function with sapply

point.count <- 50 # sets how many of the last observations to use

r.seq <- seq(1,4,by=0.001) # a sequence for different rates

plot(0,0, ylab="", xlab="", cex=0, xlim=c(0,4), ylim=c(0,1))

log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=3000),point.count)

for(i in 1:length(r.seq)) points(rep(r.seq[i],point.count),log.tree.points[,i], cex=0.2, lwd=0.2)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

system.time(log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=300),point.count))

system.time(log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=3000),point.count))

system.time(log.tree.points <- tail(sapply(r.seq, log.tree, x=0.1, n=30000),point.count))

# Martin Stoppacher                                                             #
# office@martinstoppacher.com                                                   #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#################################################################################

# Martin Stoppacher #

# office@martinstoppacher.com #

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

#################################################################################

Creating sounds out of financial data.

04/09/2015Allgemein, Code, Data Analysis, Engineering, Programming, R, Soundengineering, StatisticsMartin Stoppacher

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#

install.packages("seewave")
require("seewave")
install.packages("tuneR")
require("tuneR")

#rm(list = ls(all = TRUE)) # clear current workspace #
setwd("/Users/martinstoppacher/R Analysis/3_Index Sounds/")

library("quantmod")

getSymbols("^GSPC",from=1900)

head(GSPC)

tail(GSPC)

jpeg(filename = "SP500.jpg1975-2015.jpg", width=880,height=880,res=100)
plot(Cl(GSPC),main="S&amp;P 500 Index (closing prices)")
dev.off()

summary(Cl(GSPC))

jpeg(filename = "SP500.jpg1975-1985.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["1975/1985"],main="S&amp;P 500 Index 1975-1985 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg1985-1995.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["1985/1995"],main="S&amp;P 500 Index 1985-1995 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg1995-2005.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["1995/2005"],main="S&amp;P 500 Index 1995-2005 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg2005-2015.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["2005/2015"],main="S&amp;P 500 Index 2005-2015 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg 1975-2005 4decades-new.jpg", width=880,height=880,res=100)
plot(as.numeric(Cl(GSPC)["1975/1984"]),main="S&amp;P 500 Index 1975-1985 - 4 decades (closing prices)",ylim=c(0,2100),type="l",ylab="index values",xlab="days (10 year)")
lines(as.numeric(Cl(GSPC)["1985/1994"]),col="red")
lines(as.numeric(Cl(GSPC)["1995/2004"]),col="blue")
lines(as.numeric(Cl(GSPC)["2005/2015"]),col="green")
legend("topleft", legend = c("1975-1984","1985-1994","1995-2004","2005-2014") , lty = 1, col = c("black","red","blue","green"))
dev.off()

jpeg(filename = "SP500.jpg 1975-2005 4decades percent-new2.jpg", width=880,height=880,res=100)
plot((as.numeric(Cl(GSPC)["1975/1984"])/as.numeric(Cl(GSPC)["1975/1984"][1])-1),main="S&amp;P 500 Index 1975-1985 - 4 decades (percent changes)",ylim=c(-0.4,2.3),type="l",ylab="index values",xlab="days (10 year)")
lines((as.numeric(Cl(GSPC)["1985/1994"])/as.numeric(Cl(GSPC)["1985/1994"][1])-1),col="red")
lines((as.numeric(Cl(GSPC)["1995/2004"])/as.numeric(Cl(GSPC)["1995/2004"][1])-1),col="blue")
lines((as.numeric(Cl(GSPC)["2005/2014"])/as.numeric(Cl(GSPC)["2005/2014"][1])-1),col="green")
legend("topleft", legend = c("1975-1984","1985-1994","1995-2004","2005-2014") , lty = 1, col = c("black","red","blue","green"))
dev.off()

jpeg(filename = "SP500.jpg 1975-2005 4otherdecades percent2-new.jpg", width=880,height=880,res=100)
plot(as.numeric(Cl(GSPC)["1980/1989"])/as.numeric(Cl(GSPC)["1980/1989"][1]),main="S&amp;P 500 Index 1975-2015 - new truncation - (percent changes)",ylim=c(0.6,4.2),type="l",ylab="index values",xlab="days (10 year)",col="yellow")
lines(as.numeric(Cl(GSPC)["1975/1979"])/as.numeric(Cl(GSPC)["1975/1979"][1]),col="red")
lines(as.numeric(Cl(GSPC)["1990/1999"])/as.numeric(Cl(GSPC)["1990/1999"][1]),col="black")
lines(as.numeric(Cl(GSPC)["2000/2009"])/as.numeric(Cl(GSPC)["2000/2009"][1]),col="blue")
lines(as.numeric(Cl(GSPC)["2010/2014"])/as.numeric(Cl(GSPC)["2010/2014"][1]),col="green")
legend("topleft", legend = c("1975-1979","1980-1989","1990-1999","2000-2009","2010-2014") , lty = 1, col = c("red","yellow","black","blue","green"))
dev.off()

library("PerformanceAnalytics")

Cl(GSPC)["2010/2015"]/as.numeric(Cl(GSPC)["2005/2015"][1])
charts.PerformanceSummary(,main="",xlab="")

# percent

jpeg(filename = "SP500.jpg1975-1985-percent.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["1975/1985"]/as.numeric(Cl(GSPC)["1975/1985"][1]),main="S&amp;P 500 Index 1975-1985 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg1985-1995-percent.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["1985/1995"]/as.numeric(Cl(GSPC)["1985/1995"][1]),main="S&amp;P 500 Index 1985-1995 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg1995-2005-percent.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["1995/2005"]/as.numeric(Cl(GSPC)["1995/2005"][1]),main="S&amp;P 500 Index 1995-2005 (closing prices)")
dev.off()

jpeg(filename = "SP500.jpg2005-2015-percent.jpg", width=880,height=880,res=100)
plot(Cl(GSPC)["2005/2015"]/as.numeric(Cl(GSPC)["2005/2015"][1]),main="S&amp;P 500 Index 2005-2015 (closing prices)")
dev.off()

GSPC.cl.close &lt;- diff(Cl(GSPC))
tail(GSPC.cl.close,20)
jpeg(filename = "SP500-first-difference-example.jpg", width=880,height=880,res=100)
plot(tail(GSPC.cl.close,200),main="S&amp;P 500 first difference")
dev.off()

##

GSPC.cl.close.roc &lt;- ROC(Cl(GSPC))
tail(GSPC.cl.close.roc,5)

#prices &lt;- Cl(GSPC) # ROC is log diff!
#log_returns &lt;- diff(log(prices), lag=1)
#tail(log_returns)

jpeg(filename = "SP500-first-roc-example.jpg", width=880,height=880,res=100)
plot(tail(GSPC.cl.close.roc,200),main="S&amp;P 500 first difference")
dev.off()

##

dax.roc &lt;- na.omit(ROC(Cl(GSPC)))*100
plot(head(dax.roc,20))
plot(dax.roc)

# standard

if(abs(max(dax.roc))&gt;abs(min(dax.roc))){
dax.roc.standard &lt;- as.numeric(dax.roc/max(dax.roc))
}else{
dax.roc.standard &lt;- as.numeric(dax.roc/abs(min(dax.roc)))
}

plot(dax.roc.standard,type="l")

w&lt;-dax.roc.standard
f=41000
savewav(w,f=f ,filename = "xyz.wav")
aw&lt;-readWave("xyz.wav")
play(aw)

f=10000
savewav(w,f=f ,filename = "xyz.wav")
aw&lt;-readWave("xyz.wav")
play(aw)

f=5000
savewav(w,f=f ,filename = "xyz.wav")
aw&lt;-readWave("xyz.wav")
play(aw)

dax.roc &lt;- as.numeric(dax.roc)
dax.roc2 &lt;- NULL
for(i in 1:length(dax.roc)){
dax.roc2 &lt;- rbind(dax.roc2,((dax.roc[i]+dax.roc[i+1])/2))
}
lines &lt;- NULL
for(i in 1:length(dax.roc)){
line &lt;- rbind(dax.roc[i],dax.roc2[i])
lines &lt;- rbind(lines,line)
}
dax.roc &lt;- na.omit(lines)
tail(dax.roc)

dax.roc &lt;- na.omit(ROC(SMA(Cl(GSPC),n=500)))*100
dax.roc &lt;- na.omit(ROC(Cl(GDAXI)))*100

dax.roc.standard &lt;- as.numeric(dax.roc/max(dax.roc))
dax.roc.standard &lt;- as.numeric(dax.roc/min(dax.roc))

hist(dax.roc.standard)

w&lt;-na.omit(SMA(dax.roc.standard,n=100))
w&lt;-dax.roc.standard

for(i in 1:5){
w&lt;-c(w,w)
}

f=32000
savewav(w,f=f ,filename = "xyz.wav")
aw&lt;-readWave("xyz.wav")
play(aw)

# Martin Stoppacher #
# office@martinstoppacher.com #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#################################################################################

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

install.packages("seewave")

require("seewave")

install.packages("tuneR")

require("tuneR")

#rm(list = ls(all = TRUE)) # clear current workspace #

setwd("/Users/martinstoppacher/R Analysis/3_Index Sounds/")

library("quantmod")

getSymbols("^GSPC",from=1900)

head(GSPC)

tail(GSPC)

jpeg(filename = "SP500.jpg1975-2015.jpg", width=880,height=880,res=100)

plot(Cl(GSPC),main="S&P 500 Index (closing prices)")

dev.off()

summary(Cl(GSPC))

jpeg(filename = "SP500.jpg1975-1985.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["1975/1985"],main="S&P 500 Index 1975-1985 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg1985-1995.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["1985/1995"],main="S&P 500 Index 1985-1995 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg1995-2005.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["1995/2005"],main="S&P 500 Index 1995-2005 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg2005-2015.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["2005/2015"],main="S&P 500 Index 2005-2015 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg 1975-2005 4decades-new.jpg", width=880,height=880,res=100)

plot(as.numeric(Cl(GSPC)["1975/1984"]),main="S&P 500 Index 1975-1985 - 4 decades (closing prices)",ylim=c(0,2100),type="l",ylab="index values",xlab="days (10 year)")

lines(as.numeric(Cl(GSPC)["1985/1994"]),col="red")

lines(as.numeric(Cl(GSPC)["1995/2004"]),col="blue")

lines(as.numeric(Cl(GSPC)["2005/2015"]),col="green")

legend("topleft", legend = c("1975-1984","1985-1994","1995-2004","2005-2014") , lty = 1, col = c("black","red","blue","green"))

dev.off()

jpeg(filename = "SP500.jpg 1975-2005 4decades percent-new2.jpg", width=880,height=880,res=100)

plot((as.numeric(Cl(GSPC)["1975/1984"])/as.numeric(Cl(GSPC)["1975/1984"][1])-1),main="S&P 500 Index 1975-1985 - 4 decades (percent changes)",ylim=c(-0.4,2.3),type="l",ylab="index values",xlab="days (10 year)")

lines((as.numeric(Cl(GSPC)["1985/1994"])/as.numeric(Cl(GSPC)["1985/1994"][1])-1),col="red")

lines((as.numeric(Cl(GSPC)["1995/2004"])/as.numeric(Cl(GSPC)["1995/2004"][1])-1),col="blue")

lines((as.numeric(Cl(GSPC)["2005/2014"])/as.numeric(Cl(GSPC)["2005/2014"][1])-1),col="green")

legend("topleft", legend = c("1975-1984","1985-1994","1995-2004","2005-2014") , lty = 1, col = c("black","red","blue","green"))

dev.off()

jpeg(filename = "SP500.jpg 1975-2005 4otherdecades percent2-new.jpg", width=880,height=880,res=100)

plot(as.numeric(Cl(GSPC)["1980/1989"])/as.numeric(Cl(GSPC)["1980/1989"][1]),main="S&P 500 Index 1975-2015 - new truncation - (percent changes)",ylim=c(0.6,4.2),type="l",ylab="index values",xlab="days (10 year)",col="yellow")

lines(as.numeric(Cl(GSPC)["1975/1979"])/as.numeric(Cl(GSPC)["1975/1979"][1]),col="red")

lines(as.numeric(Cl(GSPC)["1990/1999"])/as.numeric(Cl(GSPC)["1990/1999"][1]),col="black")

lines(as.numeric(Cl(GSPC)["2000/2009"])/as.numeric(Cl(GSPC)["2000/2009"][1]),col="blue")

lines(as.numeric(Cl(GSPC)["2010/2014"])/as.numeric(Cl(GSPC)["2010/2014"][1]),col="green")

legend("topleft", legend = c("1975-1979","1980-1989","1990-1999","2000-2009","2010-2014") , lty = 1, col = c("red","yellow","black","blue","green"))

dev.off()

library("PerformanceAnalytics")

Cl(GSPC)["2010/2015"]/as.numeric(Cl(GSPC)["2005/2015"][1])

charts.PerformanceSummary(,main="",xlab="")

# percent

jpeg(filename = "SP500.jpg1975-1985-percent.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["1975/1985"]/as.numeric(Cl(GSPC)["1975/1985"][1]),main="S&P 500 Index 1975-1985 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg1985-1995-percent.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["1985/1995"]/as.numeric(Cl(GSPC)["1985/1995"][1]),main="S&P 500 Index 1985-1995 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg1995-2005-percent.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["1995/2005"]/as.numeric(Cl(GSPC)["1995/2005"][1]),main="S&P 500 Index 1995-2005 (closing prices)")

dev.off()

jpeg(filename = "SP500.jpg2005-2015-percent.jpg", width=880,height=880,res=100)

plot(Cl(GSPC)["2005/2015"]/as.numeric(Cl(GSPC)["2005/2015"][1]),main="S&P 500 Index 2005-2015 (closing prices)")

dev.off()

GSPC.cl.close <- diff(Cl(GSPC))

tail(GSPC.cl.close,20)

jpeg(filename = "SP500-first-difference-example.jpg", width=880,height=880,res=100)

plot(tail(GSPC.cl.close,200),main="S&P 500 first difference")

dev.off()

GSPC.cl.close.roc <- ROC(Cl(GSPC))

tail(GSPC.cl.close.roc,5)

#prices <- Cl(GSPC) # ROC is log diff!

#log_returns <- diff(log(prices), lag=1)

#tail(log_returns)

jpeg(filename = "SP500-first-roc-example.jpg", width=880,height=880,res=100)

plot(tail(GSPC.cl.close.roc,200),main="S&P 500 first difference")

dev.off()

dax.roc <- na.omit(ROC(Cl(GSPC)))*100

plot(head(dax.roc,20))

plot(dax.roc)

# standard

if(abs(max(dax.roc))>abs(min(dax.roc))){

dax.roc.standard <- as.numeric(dax.roc/max(dax.roc))

}else{

dax.roc.standard <- as.numeric(dax.roc/abs(min(dax.roc)))

}

plot(dax.roc.standard,type="l")

w<-dax.roc.standard

f=41000

savewav(w,f=f ,filename = "xyz.wav")

aw<-readWave("xyz.wav")

play(aw)

f=10000

savewav(w,f=f ,filename = "xyz.wav")

aw<-readWave("xyz.wav")

play(aw)

f=5000

savewav(w,f=f ,filename = "xyz.wav")

aw<-readWave("xyz.wav")

play(aw)

dax.roc <- as.numeric(dax.roc)

dax.roc2 <- NULL

for(i in 1:length(dax.roc)){

dax.roc2 <- rbind(dax.roc2,((dax.roc[i]+dax.roc[i+1])/2))

}

lines <- NULL

for(i in 1:length(dax.roc)){

line <- rbind(dax.roc[i],dax.roc2[i])

lines <- rbind(lines,line)

}

dax.roc <- na.omit(lines)

tail(dax.roc)

dax.roc <- na.omit(ROC(SMA(Cl(GSPC),n=500)))*100

dax.roc <- na.omit(ROC(Cl(GDAXI)))*100

dax.roc.standard <- as.numeric(dax.roc/max(dax.roc))

dax.roc.standard <- as.numeric(dax.roc/min(dax.roc))

hist(dax.roc.standard)

w<-na.omit(SMA(dax.roc.standard,n=100))

w<-dax.roc.standard

for(i in 1:5){

w<-c(w,w)

}

f=32000

savewav(w,f=f ,filename = "xyz.wav")

aw<-readWave("xyz.wav")

play(aw)

# Martin Stoppacher #

# office@martinstoppacher.com #

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

#################################################################################