Monday, 24 December 2018

Finds optimal alignment between two
time series x and y, and dtwDist(mx, my=mx, ...) or dist(mx, my=mx, method="DTW", ...)
calculates the distances between time series mx and my
> idx <- seq(0, 2*pi, len=100)
> a <- sin(idx) + runif(100)/10
> b <- cos(idx)

> align <- dtw(a, b, step=asymmetricP1, keep=T)

> dtwPlotTwoWay(align)


Wednesday, 19 December 2018

Gaussian white noise series

Gaussian white noise series (top) and the three-point moving average of the Gaussian white noise series (bottom)

> w = rnorm(500,0,1)
> w = rnorm(500,0,1)
> par(mfrow=c(2,1))
> plot.ts(w, main="white noise")
> v = filter(w, sides=2, rep(1/3,3))
> plot.ts(v, main="moving average")
https://www.mathclasstutor.online


Friday, 14 December 2018

One-way analysis of means

> x = c(5, 5, 5, 13, 7, 11, 11, 9, 8, 9)
> y = c(11, 8, 4, 5, 9, 5, 10, 5, 4, 10)
>
> boxplot(x,y)

> boxplot(x,y)
> y=rnorm(1000)
> f=factor(rep(1:10,100))
> boxplot(y ~ f,main="Boxplot of normal random data with model notation")

> x = rnorm(100)
> y = factor(rep(1:10,10))
> stripchart(x ~ y)
> par(mfrow=c(1,3))
> data(InsectSprays)
> boxplot(count ~ spray, data = InsectSprays, col = "lightgray")
> plot(x,y)



> x = 1:10
> y = sample(1:100,10)
> z = x+y
> lm(z ~ x+y)

Call:
lm(formula = z ~ x + y)

Coefficients:
(Intercept)            x            y 
  2.472e-14    1.000e+00    1.000e+00 

> z = x+y + rnorm(10,0,10)
> lm(z ~ x+y)

Call:
lm(formula = z ~ x + y)

Coefficients:
(Intercept)            x            y 
     -5.888        0.484        1.180 

> lm(z ~ x+y -1)

Call:
lm(formula = z ~ x + y - 1)

Coefficients:
       x         y 
-0.05044   1.13803 

> summary(lm(z ~ x+y ))

Call:
lm(formula = z ~ x + y)

Residuals:
    Min      1Q  Median      3Q     Max
-12.993  -6.778   0.519   6.582  13.599

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -5.8883     9.1775  -0.642    0.542   
x             0.4840     1.1287   0.429    0.681   
y             1.1803     0.1011  11.675 7.64e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.19 on 7 degrees of freedom
Multiple R-squared:  0.9514, Adjusted R-squared:  0.9375
F-statistic: 68.49 on 2 and 7 DF,  p-value: 2.534e-05

> x = c(4,3,4,5,2,3,4,5)
> y = c(4,4,5,5,4,5,4,4)
> z = c(3,4,2,4,5,5,4,4)
> scores = data.frame(x,y,z)
> boxplot(scores)
> scores = stack(scores)
> names(scores)
[1] "values" "ind" 
> oneway.test(values ~ ind, data=scores, var.equal=T)

One-way analysis of means

data:  values and ind
F = 1.1308, num df = 2, denom df = 21, p-value = 0.3417

> df = stack(data.frame(x,y,z))
> oneway.test(values ~ ind, data=df,var.equal=T)

One-way analysis of means

data:  values and ind
F = 1.1308, num df = 2, denom df = 21, p-value = 0.3417

> anova(lm(values ~ ind, data=df))
Analysis of Variance Table

Response: values
          Df Sum Sq Mean Sq F value Pr(>F)
ind        2   1.75 0.87500  1.1308 0.3417
Residuals 21  16.25 0.77381             
>
> kruskal.test(values ~ ind, data=df)

Kruskal-Wallis rank sum test

data:  values by ind
Kruskal-Wallis chi-squared = 1.9387, df = 2, p-value = 0.3793

Thursday, 13 December 2018

Scatter Plot with 10,000 Observations

> n <- 10000
> c1 <- matrix(rnorm(n, mean=0, sd=.5), ncol=2)
> c2 <- matrix(rnorm(n, mean=3, sd=2), ncol=2)
> mydata <- rbind(c1, c2)
> mydata <- as.data.frame(mydata)
> names(mydata) <- c("x", "y")
> (mydata,
> with(mydata,
+      plot(x, y, pch=19, main="Scatter Plot with 10,000 Observations"))
> with(mydata,
+      plot(x, y, pch=15, main="Scatter Plot with 10,000 Observations"))
> with(mydata,
+      smoothScatter(x, y, main="Scatterplot Colored by Smoothed Densities"))




Wednesday, 12 December 2018

Plot-R

 x <- rnorm(50)
> hist(x, probability = TRUE)
> lines(density(x))
> xval <- pretty(c(-3, 3), 50)
> lines(xval, dnorm(xval), col = "red")
> plot(density(rnorm(50), bw = 0.2), type = "l")
> plot(density(rnorm(50), bw = 0.6), type = "l")

Wednesday, 5 December 2018

Solving ODEs in R

Solving ODEs in R
> par(mfrow=c(1,2))
>
> plot(rnorm(100), main = "Graph 1", ylab = "Normal distribution")
> grid()

>
> legend(x = 40, y = -1, legend = "A legend")
> plot(rnorm(100), main = "Graph 2", type = "l")
> abline(v = 50)
>
> set.seed(656)
> x = c(rnorm(150, 0, 1), rnorm(150,9,1), rnorm(150,4.5,1))
>
> y = c(rnorm(150, 0, 1), rnorm(150,0,1), rnorm(150,5,1))
> XYdf = data.frame(x,y)
> plot(XYdf, pch=20)
> XY_sing = hclust(dist(XYdf), method="single")
>
> XYs3 = cutree(XY_sing,k=3)
>
> table(XYs3)
XYs3
  1   2   3
448   1   1
> XYs6 = cutree(XY_sing,k=6)
>
> table(XYs6)
XYs6
  1   2   3   4   5   6
148 150   1 149   1   1
>
> plot(XYdf, pch=20, col=XYs6)

Tuesday, 4 December 2018

R provides a number of commands for writing output

R provides a number of commands for writing output to a file
 r <- 0.11 # Annual interest rate
> term <- 10 # Forecast duration (in years)
> period <- 1/12 # Time between payments (in years)
> payments <- 100
> r <- 0.11 # Annual interest rate
> n <- floor(term/period)
> pension <- 0
> (i in 1:n) {

 for (i in 1:n) {
+   pension[i+1] <- pension[i]*(1 + r*period) + payments
+ }
> time <- (0:n)*period
> plot(time, pension)

> x <- seq(0, 5, by = 0.01)
> > y.upper <- 2*sqrt(x)
> y.lower <- -2*sqrt(x)
> y.max <- max(y.upper)
> y.upper <- 2*sqrt(x)
> y.max <- max(y.upper)
> plot(c(-2, 5), c(y.min, y.max), type = "n", xlab = "x", ylab = "y")
> plot(c(-2, 5), c(y.min, y.max), type = "n", xlab = "x", ylab = "y")
> y.min <- min(y.lower)
> plot(c(-2, 5), c(y.min, y.max), type = "n", xlab = "x", ylab = "y")
> lines(x, y.upper)
> lines(x, y.lower)
> abline(v=-1)

> points(1, 0)
> text(1, 0, "focus (1, 0)", pos=4)
> text(-1, y.min, "directrix x = -1", pos = 4)
> title("The parabola y^2 = 4*x")
> par(mfrow = c(2, 2), mar=c(5, 4, 2, 1))
> curve(x*sin(x), from = 0, to = 100, n = 1001)
> curve(x*sin(x), from = 0, to = 10, n = 1001)
> curve(x*sin(x), from = 0, to = 1, n = 1001)
> curve(x*sin(x), from = 0, to = 0.1, n = 1001)

> par(mfrow = c(1, 1))

Latin Cubes

What is Latin Cubes? The practical applications of Latin crops and related designs are factorial experiment.factorial experiments for me...