2015 사회과학원여름수학캠프 : 컴퓨터활용실습 (1 일 ) 1 변수함수의그래프확인하기 - R을활용하여 sin 등의그래프확인하기 # 1. f(x)= sin(1/x) f = function(x) sin(1/x) par(mfcol=c(3,1)) x = seq(-2,2,0.0001) plot(y~x, type='l', xlim=c(-2,2)) x = seq(-2,2,0.0001) plot(y~x, type='l', xlim=c(-0.5, 0.5)) x = seq(-2,2,0.00001) plot(y~x, type='l', xlim=c(-0.1, 0.1)) x = seq(-0.01,0.01,0.0000001) plot(y~x, type='l', xlim=c(-0.01, 0.01)) # 2. f(x) = x* sin(1/x) f = function(x) x*sin(1/x) par(mfcol=c(3,1)) x = seq(-2,2,0.0001) plot(y~x, type='l', xlim=c(-2,2)) x = seq(-0.5,0.5,0.0001) plot(y~x, type='l', xlim=c(-0.5, 0.5)) x = seq(-0.1,0.1,0.00001) plot(y~x, type='l', xlim=c(-0.1, 0.1)) # 2. f(x) = x* sin(1/x); f'(x)=2xsin(1/x)-cos(1/x) f = function(x) 2*x*sin(1/x)-cos(1/x) par(mfcol=c(3,1)) x = seq(-2,2,0.0001); plot(y~x, type='l', xlim=c(-2,2)) x = seq(-0.5,0.5,0.0001) plot(y~x, type='l', xlim=c(-0.5, 0.5)) x = seq(-0.1,0.1,0.00001) ; plot(y~x, type='l', xlim=c(-0.1, 0.1))
# 3. f(x) = x^2* sin(1/x) f = function(x) x^2*sin(1/x) par(mfcol=c(1,1)) x = seq(-2,2,0.0001) plot(y~x, type='l', xlim=c(-2,2)) x = seq(-0.5,0.5,0.0001) plot(y~x, type='l', xlim=c(-0.5, 0.5)) x = seq(-0.1,0.1,0.00001) plot(y~x, type='l', xlim=c(-0.1, 0.1)) - 그밖의다른사이트 https://www.desmos.com/calculator cos sin 함수 가 (0,1) 주위에서어떤모양을취하는지확인하자. www.wolframalpha.com plot(cos x, {x, -1,1}) http://sage.skku.edu/ plot(sin(1/x))
컴퓨터프로그램을활용한미분과적분 Wolfram alpha http://www.wolframalpha.com/ D(cos x) integrate(cos x) http://www.wolframalpha.com/calculators/integral-calculator/ http://sage.skku.edu/ diff(x^2, x) integrate(x^2,x) integrate(log(x),x) integrate(exp(-1/2*x^2), x, -infinity, infinity) R dx2x <- deriv(~ x^2, "x") ; dx2x mode(dx2x) x <- -1:2 eval(dx2x) And many more others. Maxima, Mathics, SymPy, etc.
R 로하는수치적적분 (Numerical Integration in R) 1. 1차원함수의적분 (integrate) 가 ) 함수정의 f = function(x) {... } 나 ) 적분하기 integrate(f, lower=, upper= ) f.norm1d =function(x) { 1/sqrt(2*pi)*exp(-1/2*x^2) } integrate(f.norm1d, lower=-inf, upper=inf) c3 = function(x) (x^2-2*x+3) integrate(c3, lower=0, upper=1) curve(f.norm1d, xlim=c(-6,6)) curve(c3, xlim=c(0,1)) 2. 2차원함수의적분 (cubature::adaptintegrate, R2Cuba::cuhre) 가 ) 함수정의 f=function(x) { x[1]^2 + x[2]^2 } 나 ) 적분하기 adaptintegrate(f, lowerlimit=c(, ), upperlimit=c(, )) cuhre(2, 1, f, lower=c(, ), upper=c(, )) e1 = function(x) 1/(1-x[1]*x[2]) e3 = function(x) (2-x[1]-x[2])^(-1/2) e4 = function(x) (3-x[1]-2*x[2])^(-1/2) e6 = function(x) abs(x[1]^2 + x[2]^2-0.25) e7 = function(x) abs(x[1]-x[2])^(1/2) adaptintegrate(e1, lowerlimit = c(0,0), upperlimit = c(1,1)) # pi^2/6 cuhre(2, 1, e1, lower=c(0,0), upper=c(1, 1)) vegas(2, 1, e1, lower=c(0,0), upper=c(1, 1)) emdbook::curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "contour") curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "wireframe") curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "rgl") curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "persp")
2015 사회과학원여름수학캠프 : 컴퓨터활용실습 (2일) R로하는수치미적분 CAS를활용하여방정식풀기방정식의수치적해풀기 R의최적화를활용한회귀분석 R 로하는수치적미분 (Numerical Derivation) library(numderiv) f = function(x) { a = x[1]; b = x[2]; c = x[3]; a^2+b*c-exp(c*a) } grad(f, x=c(1,1,1)) hessian(f, x=c(1,1,1))
R 로하는수치적적분 (Numerical Integration in R) 1. 1차원함수의적분 (integrate) 가 ) 함수정의 f = function(x) {... } 나 ) 적분하기 integrate(f, lower=, upper= ) f.norm1d =function(x) { 1/sqrt(2*pi)*exp(-1/2*x^2) } integrate(f.norm, lower=-inf, upper=inf) c3 = function(x) (x^2-2*x+3) integrate(c3, lower=0, upper=1) curve(f.norm1d, xlim=c(-6,6)) curve(c3, xlim=c(0,1)) 2. 2차원함수의적분 (cubature::adaptintegrate, R2Cuba::cuhre) 가 ) 함수정의 f=function(x) { x[1]^2 + x[2]^2 } 나 ) 적분하기 adaptintegrate(f, lowerlimit=c(, ), upperlimit=c(, )) cuhre(2, 1, f, lower=c(, ), upper=c(, )) e1 = function(x) 1/(1-x[1]*x[2]) e3 = function(x) (2-x[1]-x[2])^(-1/2) e4 = function(x) (3-x[1]-2*x[2])^(-1/2) e6 = function(x) abs(x[1]^2 + x[2]^2-0.25) e7 = function(x) abs(x[1]-x[2])^(1/2) adaptintegrate(e1, lowerlimit = c(0,0), upperlimit = c(1,1)) # pi^2/6 cuhre(2, 1, e1, lower=c(0,0), upper=c(1, 1)) vegas(2, 1, e1, lower=c(0,0), upper=c(1, 1)) emdbook::curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "contour") curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "wireframe") curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "rgl") curve3d(e1(c(x,y)), from=c(-1,-1), to=c(1,1), sys3d = "persp")
이변수함수그래프그리기 http://www.math.uri.edu/~bkaskosz/flashmo/graph3d/ wxmaxima http://gkerns.people.ysu.edu/maxima/maximaintro/maximaintro.pdf f(x, y) := x^2 + y^2; load(draw); draw3d(explicit(f(x,y), x, -3, 3, y, -3, 3)); Basic Algebra and Calculus in SAGE http://doc.sagemath.org/html/en/tutorial/tour_algebra.html 방정식풀기 2 차방정식 x=var('x') solve(x^2 + 2*x + 1, x) 계수가변수인 x, a, b, c=var('x a b c') solve([a*x^2+b*x+c==0], x) 2 원 1 차방정식 x, y = var('x, y') solve([x+y==4, x-y==2], x,y) 3 원방정식 var('x y z') eq1 = x+y==4 eq2 = x^2+y^2==5 eq3 = x*y*z= 3 solve([eq1, eq2, eq3], x, y, z) 수치해 (find_root, 범위 ) theta=var('theta') find_root(cos(theta)==sin(theta), 0, pi/2)
R 의최적화를활용한회귀분석 Regression : minimize Lasso Regression : minimize Ridge Regression : minimize reg=function(param) { b0=param[1]; b1=param[2]; b2=param[3] (b0+b1+5*b2-4)^2+(b0+b1*3+b2*4-5)^2+(b0+b1*5+b2*4-2)^2+(b0+b1^2+ 6*b2-3)^2 } reg.lasso=function(param, lambda=1) { b0=param[1]; b1=param[2]; b2=param[3] (b0+b1+5*b2-4)^2+(b0+b1*3+b2*4-5)^2+(b0+b1*5+b2*4-2)^2+(b0+b1^2+ 6*b2-3)^2+ lambda*(abs(b1)+abs(b2)) } reg.ridge=function(param, lambda=1) { b0=param[1]; b1=param[2]; b2=param[3] (b0+b1+5*b2-4)^2+(b0+b1*3+b2*4-5)^2+(b0+b1*5+b2*4-2)^2+(b0+b1^2+ 6*b2-3)^2+ lambda*(b1^2+b2^2) } optimx(par=c(b0=1, b1=1, b2=1), fn=reg, control=list(all.methods=t)) optimx(par=c(b0=1, b1=1, b2=1), fn=reg.lasso, control=list(all.methods=t)) optimx(par=c(b0=1, b1=1, b2=1), fn=reg.ridge, control=list(all.methods=t)) optimx(par=c(b0=0, b1=0, b2=0, lambda=1), fn=reg.ridge, control=list(all.methods=t)) optim(par=c(b0=0, b1=0, b2=0, lambda=0), fn=reg.ridge, method="bfgs") optim(par=c(b0=0, b1=0, b2=0, lambda=1), fn=reg.ridge, method="bfgs")
x1, x2 가주어졌을때 x1=c(1,3,5,2) x2=c(5,4,4,6) y=c(4,5,2,3) f=function(b) { b0=b[1]; b1=b[2]; b2=b[3]; sum((b0+b1*x1+b2*x2-4)^2) } optim(c(1,1,1), f) optim(c(1,1,1), f, method="bfgs") lm(y~x1+x2) # 위의결과와비교
R 과 LaGrange Multiplier 를활용하여최적화 비교 SAGE var('x y z p q') eq1 = 2*x-2*p-q==0 eq2 = 2*y+p-q==0 eq3 = 2*z-5*q==0 eq4 = 2*x+y-3==0 eq5 = 10*x-y+5*z==0 solve([eq1, eq2, eq3, eq4, eq5],x,y,z,p,q) R: LaGrange Multiplier 이용 eqs = function(param) { x= param[1]; y= param[2]; z=param[3]; p=param[4]; q=param[5]; c(2*x-2*p-q, 2*y+p-q, 2*z-5*q, 2*x+y-3, 10*x-y+5*z) } sol = nleqslv(c(1,1,1,1,1), eqs) 문제 maximize ln s.t ln -2*c*x+2*x+3*y-p-r*(a+b^2/x) 3*x-p 0-p b-q*(2*a)-r*x a-q*(-2*b)-2*b*log(x) -2*c+1/c-x^2-q*(2*c) x+y+z-10=0 a^2+b^2+c^2-9 a*x+b^2*log(x)-5
2015 사회과학원여름수학캠프 : 컴퓨터활용실습 (3 일 ) R 로하는행렬연산 행렬만들기 A=matrix(c(1,2,3,4), nrow=2, ncol=2, byrow=f) B=matrix(c(1,2,3,4), nrow=2, ncol=2, byrow=t) C=matrix(c(1,2,3,4), nrow=4, ncol=1, byrow=f) D=matrix(c(1,2,3,4), nrow=1, ncol=4, byrow=t) 다음행렬을변수 E, F 에저장해보자., 주의! C,D,F 는이미이름이배정되어있다. 따라서 C1, D1, E1 으로저장하자. 행렬의덧셈, 뺄셈, 곱셈, 실수배, 역행렬 A+B A-B A-C1 A %*% B C1 %*% D1 C1 %*% A 3*A solve(a) 전치행렬 t(a) 대각합 psych::tr(a) # library(psych) 항등행렬 diag(3) diag(rep(1,5)) diag(12) 대칭행렬확인 all(a == t(a)) # A=1/2*(A+t(A)) # A=t(A) %*% A; B=B %*% t(b)
정방행렬확인 (nrow(a)==ncol(a)) 대각행렬 diag(c(1,3,2,4)) diag(c(1,1,1,3,3,3)) 직교행렬확인 (all(t(q) %*% Q == diag(3))) 참고 http://www.statmethods.net/advstats/matrix.html
SAGE 1. 가모두 의하삼각행렬일때, 도하삼각행렬인가? SAGE> var('a b c d e f g h i j k l') A=matrix([[a,0,0],[b,d,0],[c,e,f]]) B=matrix([[g,0,0],[h,j,0],[i,k,l]]) A*B 2. SAGE> A=matrix(QQ,[[1,-3,0,-5],[3,-12,-2,-27],[-2,10,2,24],[-1,6,1,14]]); b=matrix(qq, 4, 1, [-7,-33,29,17]); C=A.augment(b); C.echelon_form() 3. SAGE> A=matrix(QQ,[[0,0,0,1,2,-1],[1,2,0,0,1,-1],[1,2,2,0,-1,1]]); b=matrix(qq, 3, 1, [2,0,2]); C=A.augment(b); C.echelon_form() 4. SAGE> A=matrix(QQ, [[1,3,1,6],[2,2,3,3],[-1,1,0,1]]); b=matrix(qq, 3, 1, [8,14,-2]); C=A.augment(b); C.echelon_form()
5. SAGE> A=matrix(QQ, [[1,2,2,8,0,2,2],[2,4,1,7,1,-1,-5],[-1,-2,1,1,1,0,3],[3,6,3,15,0,3,0]]) b=matrix(qq, 4, 1, [3,0,-2,6]) C=A.augment(b); C.echelon_form()
2015 사회과학원여름수학캠프 : 컴퓨터활용실습 (4 일 ) Visualizing 2-variable function #http://stackoverflow.com/questions/20549540/how-to-create-3d-matlab-style-s urface-plots-in-r for (r_ in c(1,0, 0.5, 0.2, 0.1)) { points = seq(-r_-0.1, r_-0.1, length=20) #create a grid XY = expand.grid(x=points,y=-points) # A z-function Zf <- f # populate a surface Z <- Zf(XY$X, XY$Y) zlim <- range(z) zlen <- zlim[2] - zlim[1] + 1 jet.colors <- # function from grdevices package colorramppalette(c("#00007f", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) colorzjet <- jet.colors(100) # 100 separate color require(rgl) open3d() rgl.surface(x=points, y=matrix(z,20), coords=c(1,3,2),z=-points, color=colorzjet[ findinterval(z, seq(min(z), max(z), length=100))] ) axes3d() rgl.snapshot("copymatlabstyle.png") }
covariance matrix of linear transfomation of random vectors X <- matrix(rnorm(3*n.subj), nrow=3) A <- matrix(c(0.5, 0.4, -0.3, 0.1, -0.6, 0.1, 0.5, 0.1, -0.2, 0.05, -0.5, -0.1, 0.5, 0.1, 0.4, 0.4, -0.15, 0.01, 0.1, 0.0, 0.1), ncol=3, byrow=t) sigma1 = cov(t(a %*% X)) sigma2 = A %*% cov(t(x)) %*% t(a) sigma1 == sigma2 all.equal(sigma1, sigma2) # Almost true?
PCA n.subj <- 100 pc <- matrix(rnorm(3*n.subj), nrow=3) A <- matrix(c(0.5, 0.4, -0.3, 0.1, -0.6, 0.1, 0.5, 0.1, -0.2, 0.05, -0.5, -0.1, 0.5, 0.1, 0.4, 0.4, -0.15, 0.01, 0.1, 0.0, 0.1), ncol=3, byrow=t) dat <- t(a %*% pc) df.dat <- data.frame(dat) plot(df.dat) dat.prcomp <- prcomp(dat, cor=f ) par(mfcol=c(3,3)) cor(pc[1,],dat.prcomp$x[,"pc1"]); plot(pc[1,],dat.prcomp$x[,"pc1"]) cor(pc[2,],dat.prcomp$x[,"pc1"]); plot(pc[2,],dat.prcomp$x[,"pc1"]) cor(pc[2,],dat.prcomp$x[,"pc1"]); plot(pc[3,],dat.prcomp$x[,"pc1"]); cor(pc[1,],dat.prcomp$x[,"pc2"]) cor(pc[2,],dat.prcomp$x[,"pc2"]) cor(pc[3,],dat.prcomp$x[,"pc2"]) cor(pc[1,],dat.prcomp$x[,"pc3"]) cor(pc[2,],dat.prcomp$x[,"pc3"]) cor(pc[3,],dat.prcomp$x[,"pc3"]) factanal(dat, factors=3) solve(cor(dat)) det(cor(dat))
eigen decomposition A.eigen <- eigen(a) A.eigen$val A.eigen$vec # 대칭행렬 A 가양의정부호 (positive definite) 행렬임을확인하는방법은? # 교재 45 쪽의 행렬을구성해보자. # 은대칭행렬인가? # 은정사영행렬인가?
PCA 와 MDS 비교 # from http://www.statmethods.net/advstats/mds.html library(mass) S=matrix(c(1,0.7,0.7,0.8), ncol=2) dat <- mvrnorm(n=1000, mu=c(0,0), S) d <- dist(dat) # euclidean distances between the rows fit <- cmdscale(d,eig=true, k=1) # k is the number of dim fit # view results MDS1 <- fit$points[,1] MDS2 <- fit$points[,2] dat.scale <- scale(dat, scale=f) apply(dat.scale, 2, mean) cov.dat.scale <- t(dat.scale) %*% dat.scale eigen(cov.dat.scale) dat.pca <- prcomp(dat, center = TRUE, scale. = F) PC1 <- dat.pca$x[,"pc1"] PC2 <- dat.pca$x[,"pc2"] plot(mds1, PC1) plot(dat[,1], dat[,2]) plot(mds1, MDS2) plot(pc1, PC2)