Q1
#This is a deom of brownian motion, which is the random motion of atoms or molecules in the gas or liquid.
#It will plot a lot pictures so I close the ouput and represent code only.
brownian_Motion <- function(n = 11, pause = 0.05, nIter = 100, ...) {
x = rnorm(n)
y = rnorm(n)
i = 1
repeat {
plot(x, y, ...)
text(x, y, cex=0.5)
x = x + rnorm(n)
y = y + rnorm(n)
Sys.sleep(pause)
i = i + 1
if(i == nIter)
break
}
}
###
## test it
brownian_Motion(xlim = c(-20, 20), ylim = c(-20, 20),
pch = 21, cex = 2, col = "cyan", bg = "lavender")
Q2
# median minimizes sum of absolute deviation
#
cal_median <- function(x) {
n <- length(x)
h <- round((n+1)/2)
m <- unique(apply(expand.grid(x[1:h], x[h:n]), 1, sum)/2)
sad <- rep(NA, length(m))
sad <- sapply(1:length(m),function(i){sad[i] = sum(abs(x - m[i]))})
ll <- length(which(sad == min(sad)))
if( ll == 1)
m[which(sad == min(sad))]
else
sort(m[which(sad == min(sad))])[(ll+1)/2]
}
###
##
# test it
#
cal_median(women$height)
## [1] 65
###
Q3
dta3 <- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/dataM/Data/hs0.txt",header = T)
#filter the students whose race are asian
dta.asian <- subset(dta3, race=="asian")
#find the pearson correlation
r0 <- cor(dta.asian$math, dta.asian$socst)
nIter <- 1001
cnt <- rep(0,nIter)
#permutation
sum(sapply(1:nIter, function(i){new <- sample(dta.asian$read)
r <- cor(new, dta.asian$math)
ifelse(r0 <= r,cnt[i] <- 1,cnt[i] <- 0)
}))/nIter
## [1] 0.03496503
Q4
# riffle
# interlace two lists of numbers, starting with the longer list,
# without repeating
#
riffle <- function(x, y) {
if (length(x) >= length(y)){
a <- x
b <- y
}else{
a <- y
b <- x
}
dl <- length(a) - length(b)
b[(length(b)+1):(length(b)+dl)] <- rep(NA, dl)
if (length(x) >= length(y)){
z <- c(na.omit(as.numeric(t(cbind(a,b)))))
}else{
z <- c(na.omit(as.numeric(t(cbind(b,a)))))
}
return(z)
}
###
##
# test it
riffle(1:2, 50:55)
## [1] 1 50 2 51 52 53 54 55
###
Q5
#
nameCombo <- function( x ) {
n <-length(x)
#coerce argument to a list
ll <- as.list(x)
for(i in 2:length(x)) {
#find all combination of the elements of x taken i at a time
indices <- combn(1:n, i)
#combine the original list and new combination
for(j in 1:dim(indices)[2]) {
new_ll <- list(x[indices[, j]])
ll <- c(ll, new_ll)
}
}
return(ll)
}
#
#test it
nameCombo(c("V1", "V2", "V3"))
## [[1]]
## [1] "V1"
##
## [[2]]
## [1] "V2"
##
## [[3]]
## [1] "V3"
##
## [[4]]
## [1] "V1" "V2"
##
## [[5]]
## [1] "V1" "V3"
##
## [[6]]
## [1] "V2" "V3"
##
## [[7]]
## [1] "V1" "V2" "V3"
#
Q6
#
# plot a color strip
#
plot.new()
plot.window(xlim = c(0, 6), ylim = c(0, 6), asp = 1)
my_cl <- c("indianred", "orange", "yellow", "forestgreen", "dodgerblue",
"violet", "purple")
#
for(i in 1:6){
n = 7-i
m = -2+i
for(j in 1:n){
rect(j+m, j-1, j+m+1, j, col = my_cl[i])
}
}
#
for(i in 1:5){
col = i+2
m = 5-i
for(j in 1:i){
rect(j-1, j+m, j, j+m+1, col = my_cl[col])
}
}
Q7
cointoss <- function(P = c(0.5,0.5)){
x <- sample(c("H","T"),100,replace = T,prob = P)
length <- max(rle(x)[[1]])
value <- rle(x)[[2]][which(rle(x)[[1]] == max(rle(x)[[1]]))]
return(c(length,value))
}
vec <- seq(0,1,0.1)
nvec <- 1-vec
for(i in 1:length(vec)){
tmp <- cointoss(P = c(vec[i],nvec[i]))
print(tmp)
}
## [1] "100" "T"
## [1] "21" "T"
## [1] "12" "T"
## [1] "14" "T"
## [1] "7" "T"
## [1] "9" "T"
## [1] "8" "H"
## [1] "10" "H"
## [1] "10" "H"
## [1] "27" "H"
## [1] "100" "H"
Q8
dta8 <-scan("http://titan.ccunix.ccu.edu.tw/~psycfs/dataM/Data/plasma.txt",what = "character",sep = " ")
is.na(dta8) <- dta8[1]
dta8L <- rep(NA,length(dta8))
for(i in 1:length(dta8)){
if(dta8[i] == dta8[1]){
dta8L[i] <- NA
}else{
dta8L[i] = dta8[i]
}
}
## Error in if (dta8[i] == dta8[1]) {: 需要 TRUE/FALSE 值的地方有缺值
dta8 <- na.omit(dta8L)
dta8 <- as.numeric(dta8)
for(i in 1:length(dta8)){
if (dta8[i] < 0 ){
dta8[i] <- NA
}else{
dta8[i] <- dta8[i]
}
}
dta8 <- as.data.frame(matrix(dta8,84,9))
colnames(dta8) <- c("plasma_ascorbic_acid","leucocyte_ascorbic_acid",
"whole_blood_ascorbic_acid","thiamin_status",
"grip_strength","red_cell_transketolase",
"reaction_time","folate_serum","folate_red_cell")
str(dta8)
## 'data.frame': 84 obs. of 9 variables:
## $ plasma_ascorbic_acid : num 22 0 103 67 75 65 59 46 16 37 ...
## $ leucocyte_ascorbic_acid : num 48 32 0 0 0 21 4 14 11 18 ...
## $ whole_blood_ascorbic_acid: num 56 126 55 71 54 64 60 57 58 176 ...
## $ thiamin_status : num 54 42 133 132 130 74 56 43 58 57 ...
## $ grip_strength : num 110 45 32 0 9 10 0 9 10 12 ...
## $ red_cell_transketolase : num 50 55 40 58 38 49 36 60 52 NA ...
## $ reaction_time : num 31 54 124 56 77 28 40 63 70 71 ...
## $ folate_serum : num 5 49 12 0 11 24 13 35 30 27 ...
## $ folate_red_cell : num 180 70 80 60 70 NA 100 56 32 93 ...
for(i in 1:12){
if (i <= 9){
dta8$patient[(1+(i-1)*7):(i*7)] <- rep(paste0("S",10,1),7)
}else{
dta8$patient[(1+(i-1)*7):(i*7)] <- rep(paste0("S",1,i),7)
}
}
dta8$week <- rep(c(1,2,6,10,14,15,16),12)
dta8 <- dta8[,c(10,11,1,2,3,4,5,6,7,8,9)]
str(dta8)
## 'data.frame': 84 obs. of 11 variables:
## $ patient : chr "S101" "S101" "S101" "S101" ...
## $ week : num 1 2 6 10 14 15 16 1 2 6 ...
## $ plasma_ascorbic_acid : num 22 0 103 67 75 65 59 46 16 37 ...
## $ leucocyte_ascorbic_acid : num 48 32 0 0 0 21 4 14 11 18 ...
## $ whole_blood_ascorbic_acid: num 56 126 55 71 54 64 60 57 58 176 ...
## $ thiamin_status : num 54 42 133 132 130 74 56 43 58 57 ...
## $ grip_strength : num 110 45 32 0 9 10 0 9 10 12 ...
## $ red_cell_transketolase : num 50 55 40 58 38 49 36 60 52 NA ...
## $ reaction_time : num 31 54 124 56 77 28 40 63 70 71 ...
## $ folate_serum : num 5 49 12 0 11 24 13 35 30 27 ...
## $ folate_red_cell : num 180 70 80 60 70 NA 100 56 32 93 ...