Functions_Homework
Li-Ping Yen
2018.4.9
Q1
NT_CVT_US<- function() {
M<- readline(prompt="NT:")
cat("US:$",as.numeric(M)/29.32,"", sep = "")
}
NT_CVT_US()## NT:
## US:$NAQ2
library(datasets)
dta2<-ChickWeight
head(dta2)## Grouped Data: weight ~ Time | Chick
## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
## 6 93 10 1 1str(dta2)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 578 obs. of 4 variables:
## $ weight: num 42 51 59 64 76 93 106 125 149 171 ...
## $ Time : num 0 2 4 6 8 10 12 14 16 18 ...
## $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Diet : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "outer")=Class 'formula' language ~Diet
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "labels")=List of 2
## ..$ x: chr "Time"
## ..$ y: chr "Body weight"
## - attr(*, "units")=List of 2
## ..$ x: chr "(days)"
## ..$ y: chr "(gm)"sapply(split(dta2,dta2$Chick),
function(x) lm(weight ~Time, data = x)$coef)## 18 16 15 13 9 20 10
## (Intercept) 39 43.392857 46.83333 43.384359 52.094086 37.667826 38.695054
## Time -2 1.053571 1.89881 2.239601 2.663137 3.732718 4.066102
## 8 17 19 4 6 11
## (Intercept) 43.727273 43.030706 31.21222 32.86568 44.123431 47.921948
## Time 4.827273 4.531538 5.08743 6.08864 6.378006 7.510967
## 3 1 12 2 5 14
## (Intercept) 23.17955 24.465436 21.939797 24.724853 16.89563 20.52488
## Time 8.48737 7.987899 8.440629 8.719861 10.05536 11.98245
## 7 24 30 22 23 27
## (Intercept) 5.842535 53.067766 39.109666 40.082590 38.428074 29.858569
## Time 13.205264 1.207533 5.898351 5.877931 6.685978 7.379368
## 28 26 25 29 21 33
## (Intercept) 23.984874 20.70715 19.65119 5.882771 15.56330 45.830283
## Time 9.703676 10.10316 11.30676 12.453487 15.47512 5.855241
## 37 36 31 39 38 32
## (Intercept) 29.608834 25.85403 19.13099 17.03661 10.67282 13.69173
## Time 6.677053 9.99047 10.02617 10.73710 12.06051 13.18091
## 40 34 35 44 45 43
## (Intercept) 10.83830 5.081682 4.757979 44.909091 35.673121 52.185751
## Time 13.44229 15.000151 17.258811 6.354545 7.686432 8.318863
## 41 47 49 46 50 42
## (Intercept) 39.337922 36.489790 31.662986 27.771744 23.78218 19.86507
## Time 8.159885 8.374981 9.717894 9.738466 11.33293 11.83679
## 48
## (Intercept) 7.947663
## Time 13.714718Q3
L<- function(df){
curve(dnorm(x),-4,4,col= 2, ylab = "dnorm(x)", lwd = 2)
for(i in 1:length(df))curve(dt(x,df[[i]]),col=3,lty=2,add=TRUE)
}
L(df=1:32)
Q4
library(pacman)
pacman::p_load(MASS, tidyverse)
# method 1 #做法簡單
aggregate( . ~ Type, data = Cushings, mean)## Type Tetrahydrocortisone Pregnanetriol
## 1 a 2.966667 2.44
## 2 b 8.180000 1.12
## 3 c 19.720000 5.50
## 4 u 14.016667 1.20# method 2 # wide format
sapply(split(Cushings[,-3], Cushings$Type), function(x) apply(x, 2, mean))## a b c u
## Tetrahydrocortisone 2.966667 8.18 19.72 14.01667
## Pregnanetriol 2.440000 1.12 5.50 1.20000# method 3 #long format
do.call("rbind", as.list(
by(Cushings, list(Cushings$Type), function(x) {
y <- subset(x, select = -Type)
apply(y, 2, mean)
}
)))## Tetrahydrocortisone Pregnanetriol
## a 2.966667 2.44
## b 8.180000 1.12
## c 19.720000 5.50
## u 14.016667 1.20# method 4 #使用dplyr裡面的pipe功能，快速分組運算
Cushings %>%
group_by(Type) %>%
summarize( t_m = mean(Tetrahydrocortisone), p_m = mean(Pregnanetriol))## # A tibble: 4 x 3
## Type t_m p_m
## <fct> <dbl> <dbl>
## 1 a 2.97 2.44
## 2 b 8.18 1.12
## 3 c 19.7 5.50
## 4 u 14.0 1.20# method 5 #查看資料屬性
Cushings %>%
nest(-Type) %>%
mutate(avg = map(data, ~ apply(., 2, mean)),
res_1 = map_dbl(avg, "Tetrahydrocortisone"),
res_2 = map_dbl(avg, "Pregnanetriol")) ## Warning: package 'bindrcpp' was built under R version 3.4.4## # A tibble: 4 x 5
## Type data avg res_1 res_2
## <fct> <list> <list> <dbl> <dbl>
## 1 a <data.frame [6 × 2]> <dbl [2]> 2.97 2.44
## 2 b <data.frame [10 × 2]> <dbl [2]> 8.18 1.12
## 3 c <data.frame [5 × 2]> <dbl [2]> 19.7 5.50
## 4 u <data.frame [6 × 2]> <dbl [2]> 14.0 1.20Q5
K<- function(n,mu,sigm){
set.seed(0221)
sample<-rnorm(n,mu,sigm)
plot(x = 1:n, y = cumsum(sample)/1:n, type = "l", col = 3,
xlab = "Sample Size", ylab = "Running Average")
abline(h = mu, col = 2, lty = 2)
}
#EX
K(3000,20,10)
Q6
setwd("/Users/tayloryen/Desktop/大學/成大課業/大四下/資料管理/0409/HW")
dta6<-read.table("q6.txt", header=T)
# function
c.stat <- function(data, n = length(data)){
cden <- 1-(sum(diff(dta6[1:n,1])^2)/(2*(n-1)*var(dta6[1:n,1])))
sc <- sqrt((n-2)/((n-1)*(n+1)))
pval <- 1-pnorm(cden/sc)
return(list(c=cden,z=cden/sc,pvalue=pval))
}
c.stat(dta6,36) ## $c
## [1] 0.660105
##
## $z
## [1] 4.073882
##
## $pvalue
## [1] 2.311794e-05Q7
library(rgl)## Warning: package 'rgl' was built under R version 3.4.4knitr::knit_hooks$set(webgl = hook_webgl)
dta <- read.table("q7.txt", header = TRUE)
ssq <- function(mu, sigma, y) {sum(((y - mu) / sigma)^2)}
vssq <- Vectorize(ssq, c("mu", "sigma"))
RR<- function(data = NA){
x <- seq(mean(data) - 5, mean(data) + 5, by = 0.1)
y <- sd(data)
z <- vssq(x, sd(data), data)
rgl::plot3d(x, y, z,
xlab = "mu", ylab = "sigma(control)", zlab = "ssq")
}
##The End
RR(data=dta$math)You must enable Javascript to view this page properly.