Skip to content

Instantly share code, notes, and snippets.

@syou6162
Created September 13, 2008 06:46
Show Gist options
  • Save syou6162/10566 to your computer and use it in GitHub Desktop.
Save syou6162/10566 to your computer and use it in GitHub Desktop.
Tsukuba.Rの資料用
library(scatterplot3d)
steepest_descent <- function(epsilon,tau,x_bar1=-10,x_bar2=-20,beta=1,n=2000){
#関数内関数の定義
f <- function(x){x1 <- x[1];x2 <- x[2];return(100*(x1-x2)^2+(x1-1)^2)}
#数式微分
trig.exp <- expression(100*(x1-x2)^2+(x1-1)^2)
D.sc1 <- D(trig.exp, "x1")
D.sc2 <- D(trig.exp, "x2")
#最急降下法の計算
x_vec <- matrix(0,n+1,2,byrow=TRUE)
for(i in seq(n)){
nabla_f <- c((function(x1,x2){eval(D.sc1)})(x_bar1,x_bar2),
(function(x1,x2){eval(D.sc2)})(x_bar1,x_bar2))
while(f(c(x_bar1,x_bar2) - beta * nabla_f) >
f(c(x_bar1,x_bar2)) - epsilon * beta * sum(c((function(x1,x2){eval(D.sc1)})(x_bar1,x_bar2),(function(x1,x2){eval(D.sc2)})(x_bar1,x_bar2))^2)){
beta <- tau * beta
}
a <- beta
x_bar1_old <- x_bar1
x_bar1 <- x_bar1 - a * (function(x1,x2){eval(D.sc1)})(x_bar1,x_bar2)
x_bar2 <- x_bar2 - a * (function(x1,x2){eval(D.sc2)})(x_bar1_old,x_bar2)
x_vec[i,] <- c(x_bar1,x_bar2)
}
#3次元のplot関係
x <- seq(-1,1,length=50)
y <- x
z <- outer(x,y,function(x,y){mapply(function(x,y){f(c(x,y))},x,y)})
# title(main="最急降下法",sub=paste(deparse(substitute(tau)),"=",tau,sep=""))
cat(deparse(substitute(x)))
persp(x,y,z,theta=150,phi=10,expand=0.5) -> res
x <- x_vec[,1]
y <- x_vec[,2]
z <- mapply(function(x,y){f(c(x,y))},x,y)
title(main="最急降下法",sub=paste(
deparse(substitute(tau)),"=",tau,
" ",
deparse(substitute(epsilon)),"=",epsilon,
sep=""))
lines(trans3d(x,y,z, pmat = res),col="red",lwd=5)
#結果を明示的な形で返さないようにする
return(invisible(x_vec))
}
plot
x_vec
plot(1:10)
tau <- 0.1
substitute(tau)
deparse(substitute(tau))
png(paste(substitute(tau),"hoge.png",sep=""))
dmp <- function(x) cat(deparse(substitute(x)),"=",x,"\n")
x <- 3
dmp(x)
x = 3
y <- "abc"
dmp(y)
dev.off()
deparse(substitute(c))
setwd("~/Desktop")
par(ask=FALSE)
mat <- matrix(1:9,3,3,byrow=TRUE)
mat
layout(mat)
for(epsilon in c(0.1,0.5,0.9)){
for(tau in c(0.1,0.5,0.9)){
steepest_descent(epsilon,tau)
}
}
getwd()
#直積を書き出す関数
write.direct.product <- function(x,y){
e <- expand.grid(x,y)
m <- mapply(function(x,y){paste("(",x,",",y,")",sep="")},e[,1],e[,2])
return(paste("{",paste(m,collapse=","),"}",sep=""))
}
write.direct.product(1:3,4:6)
#R \times RからRへの写像の例
apply(expand.grid(1:3,4:6),1,function(x){x[1]*x[2]})
apply(expand.grid(letters[1:3],1:3),1,function(d){Reduce(function(x,init){paste(init,x,sep="")},rep(d[1],d[2]),init="")})
#和集合を作る操作
my.union <- function(list){
Reduce(function(x,init){unlist(union(x,init))},list,list())
}
my.union(list(1:30,2:50,4:1,list(1,2,list(3,4))))
my.union(list(list(1,2),list(3,4)))
#共通集合を作る操作
my.intersect <- function(list){
Reduce(function(x,init){unlist(intersect(x,init))},list,list[[1]])
}
my.intersect(list(1:30,2:50,4:10))
#冪集合を返す関数。空集合はNAにしてみた
powerset <- function(list){
c(NA,unlist(lapply(seq(length(list)),function(x){apply(combn(list,x),2,function(y){y})}),recursive=FALSE))
}
#冪集合の例
powerset(list("a","b","c"))
lapply(powerset(list("a","b","c")),unlist)
powerset(list("a","b","c","d",1:10))
lapply(powerset(list("a","b","c","d")),unlist)[[x]]
my.union(mapply(function(x){lapply(powerset(list("a","b","c","d")),unlist)[[x]]},c(16,15)))
my.intersect(mapply(function(x){lapply(powerset(list("a","b","c","d")),unlist)[[x]]},c(16,1)))
write.list <- function(list){
write.vector <- function(vector){
return(paste("{",paste(vector,collapse=","),"}",sep=""))
}
write.list.intern <- function(list){
ifelse(is.list(list),Reduce(function(init,x){paste(init,write.list.intern(x),sep="")},list,""),write.vector(list))
}
gsub("NA","",paste("{",gsub("\\}\\{","\\},\\{",write.list.intern(list)),"}",sep=""))
}
write.list(list("a","b","c","d",list(1:10,list(1,2))))
write.list(lapply(powerset(list("a","b","c","d")),unlist))
a <- "Misho"
b <-"キラッ"
paste(a,b,sep="☆")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment