おまけ：日本語を表示させたいとき
課題１：素数を数える
課題２：平方根を与える組
課題３：積分の近似
課題４：級数の和
課題５：正多角形
課題６：素因数分解
課題７：階乗
課題８：円の面積の近似値
課題９：４つの円が重なる面積
課題１０：対数グラフ
課題１１：F(x)=1+1/2^x+1/3^x +...+1/100^xのグラフ
課題１２：半球
課題１３：パッケージbiOpsを使った画像処理
課題１４：ランダムな位置と色で円を描く
課題１５：パッケージrgl、および関数spheres3dについて
おまけ：日本語を表示させたいとき †
グラフを描くとき、HiraMaru?フォントなら日本語も表示できる（Macを使うとき）
quartzFonts(HiraMaru=quartzFont(rep("HiraMaruProN-W4", 4))) # Rを起動時に実行
par(family="HiraMaru") # グラフを描画する度に実行
↑
課題１：素数を数える †
素数を判定する関数
hantei=function(n){
hantei=1
if (n==2 | n==3) {
return(hantei) # 2と3は素数なので判定を1のまま終了
}
m=ceiling(sqrt(n))
for (i in 2:m){
if (n %% i==0) {
hantei=0
return(hantei) # 素数じゃないので判定を0にして終了
}}
return(hantei)
}
n 以下の素数を求める関数
prime=function(n){
z=0
prime_array=c(NA)
for (i in 2:n){
if (hantei(i)==1) { # ここでhantei関数を使う
z=z+1
prime_array[z]=i
}}
return(prime_array)
}
結果のグラフを描画する
x=c(1000, c(1:5)*10000) # xに1000,100000,200000,...,500000 を格納する
G=c(NA)
H=c(NA)
for (i in 1:length(x)){
n=x[i]
G[i]=length(prime(n)) # 計算値
H[i]=n / log(n) # 理論値
}
# 相対頻度をグラフにするため、GとHを個数xで割っている
plot(x,G/x,type="b",ylim=c(0,0.2),xlim=c(0,50000),col="blue",ylab="")
par(new="T")
plot(x,H/x,type="b",ylim=c(0,0.2),xlim=c(0,50000),col="red",lty=2,main="x以下の素数の相対頻度と理論近似値",ylab="相対頻度とその理論近似値")
legend(1000, 0.20, c("相対頻度", "理論近似値"), col = c("blue","red"), lty = c(1,2), pch = c(1,1))
↑
課題２：平方根を与える組 †
n=100
z=0
for (i in 1:n){
for (j in 1:n){
m=i^2+j^2
if (ceiling(sqrt(m))^2==m) {z=z+1;cat(i,"^2","+",j,"^2=",m,"\n")}
# ceilingは、引数の数の最も大きな整数、sqrtは平方根を求める
}}
cat("1=<x,y=<",n,"のとき平方数を与える組(x,y)の個数は",z,"\n")
↑
課題３：積分の近似 †
curve(sin(x),xlim=c(0,pi),ylim=c(0,2),col="red") # sin(x)を描画する：赤
par(new=T) # 重ね描き
n=100 # nは分割数 : 100等分する
x=seq(0,pi,length=n)
for(i in x){
polygon(c(i,i,i+pi/n,i+pi/n),c(0,sin(i),sin(i),0)) # 分割された面積を描画する：黒
}
S= sum(sin(x))*(pi/n) # 面積
cat(S) # ２に値が近いことを確認する
↑
課題４：級数の和 †
n=20
S=0 # 初期値
for( i in 1:n ){
S=S+ i^3
}
cat(S)
↑
課題５：正多角形 †
n=seq(10, 100, by=10) # 正多角形の辺の数
S=0 # 初期値
for( i in 1:length(n)){
S[i]= sin(pi/n[i])*cos(pi/n[i])*n[i]
}
plot(S, type="b", col="blue", main="正多角形の面積")
abline(h=pi,col="red",lwd=2) # 3.14のラインを引き、半径１の円の面積に近づくことを確認する
↑
課題６：素因数分解 †
課題２のhantei,prime関数を使うので、Rに読み込ませておくこと。
素因数分解する。
factor=function(n){
F=prime(n)
L=length(F)
G=c(NA)
z=0
# n以下の素数の中でｎを割り切る素数をGに格納する
for (i in 1:L){
if (n %% F[i]==0) {
z=z+1
G[z]=F[i]
}}
LL=length(G) # 素数の個数
bekiv=c(NA) # べき数
# nをわりきる素数が何冪まで割れるかを調べ、bekivに格納する
for (i in 1:LL){
m=G[i]
beki=1
N=n
while (N%%m==0){
N=N/m
beki=beki+1
}
bekiv[i]=beki-1
}
# 数式に表す
express=paste(n,"=",sep="")
for (i in 1:LL){
express=paste(express,G[i],"^",bekiv[i],sep="")
if (i<LL) express=paste(express,"*",sep="")
}
return(express)
}
使用例
# 例１
factor(2^3*3^2)
# 例２
factor(189)
↑
課題７：階乗 †
n=10
P=1 # 初期値
for(i in 1:n){
P=P*i
}
cat(P) # 結果
↑
課題８：円の面積の近似値 †
Sn=c(NA) # 初期値
for (i in 1:10){
n=i*1000
x=runif(n,-1,1)
y=runif(n,-1,1)
# 条件：which(x^2+y^2<1) を満たす点の個数を求める
Sn[i]=4*length(which(x^2+y^2<1))/n # 円の面積の近似値
}
x=c(1:10)*1000
plot(x,Sn,col="blue",type="b",main="πの乱数による推定値",lwd=2)
abline(h=pi,col="red",lwd=2)
↑
課題９：４つの円が重なる面積 †
4つの円を描画する
par(mfrow=c(2,1))
t=seq(0,2*pi,len=10000)
plot(cos(t),sin(t),xlim=c(-1,2),ylim=c(-1,2),col="blue", type="l", xlab="", ylab="", main="円が重なった領域の面積の近似値を求める")
par(new=T)
plot(cos(t),sin(t)+1,xlim=c(-1,2),ylim=c(-1,2),col="red", type="l" , xlab="", ylab="")
par(new=T)
plot(cos(t)+1,sin(t),xlim=c(-1,2),ylim=c(-1,2),col="brown", type="l", xlab="", ylab="")
par(new=T)
plot(cos(t)+1,sin(t)+1,xlim=c(-1,2),ylim=c(-1,2),col="green", type="l", xlab="", ylab="")
円が重なりあう部分についての処理
n=seq(100,1000,by=100)
Sn = rep(0, length(n))
for( i in 1:length(n)){
# 乱数を発生させる
x = runif(n[i])
y = runif(n[i])
S1=which(x^2+y^2<1)
S2=which((x-1)^2+y^2<1)
S3=which(x^2+(y-1)^2<1)
S4=which((x-1)^2+(y-1)^2<1)
# 重なった部分をSとする
S=intersect(intersect(S1,S2),intersect(S3,S4)) # intersectは積集合を求める
par(new="T")
plot(x[S],y[S],col="grey", xlim=c(-1,2),ylim=c(-1,2)) # Sを灰色に塗る
Sn[i] = length(S)/n[i] # 面積の近似を求める
cat(n[i],"個のときの面積の近似値は、", Sn[i], "です\n")
}
par(new="F")
plot(n,Sn,type="l", main="面積の近似値の折れ線グラフ")
↑
課題１０：対数グラフ †
curve(log(x),from=0,to=101,col="blue",main="log(x)のグラフ")
↑
課題１１：F(x)=1+1/2^x+1/3^x +...+1/100^xのグラフ †
F_function = function(x){
n=100
F=0 # 初期値
for(i in 1:n){
F = F + x^i
}
return(F)
}
curve(F_function,from=0, to=1,col="blue", main="F(x)=1+1/2^x+1/3^x+...+1/100^xのグラフ")
↑
課題１２：半球 †
library(scatterplot3d)
epsilon = 0
s=seq(0,pi/2,len=100)
t=seq(epsilon,2*pi,len=100)
x=as.vector(sin(s)%*%t(cos(t)))
y=as.vector(sin(s)%*%t(sin(t)))
z=as.vector(cos(s)%*%t(rep(1,100)))
scatterplot3d(x,y,z,color=rgb(.5, .8, 1, .2))
↑
課題１３：パッケージbiOpsを使った画像処理 †
パッケージbiOpsはWindows版でしか使えない。テストでは大まかな処理を理解すればOK。
library(biOps)
# 画像データをXに入れる
# X=readJpeg("/auto_mnt/home3/share/students/sakata/flower/flower.jpg") # ここ画像の格納場所によって任意
par(mfrow=c(2,2)) # 画像を2×2で配置する
# 赤の画像を作る
# 画像データを入れる配列を作る。dim(X)[1],dim(X)[2]は画像xの縦横
XR=array(dim=c(dim(X)[1],dim(X)[2],3))
XR[,,1]=imgRedBand(X) # 赤色のデータ
XR[,,2]=0 # それ以外は0にする
XR[,,3]=0
plot(imagedata(XR)) # 描画
# 緑の画像
XG=array(dim=c(dim(X)[1],dim(X)[2],3))
XG[,,1]=0
XG[,,2]=imgGreenBand(X)
XG[,,3]=0
plot(imagedata(XG))
# 青の画像
XB=array(dim=c(dim(X)[1],dim(X)[2],3))
XB[,,1]=0
XB[,,2]=0
XB[,,3]=imgBlueBand(X)
plot(imagedata(XB))
plot(imagedata(X)) # 元の画像
↑
課題１４：ランダムな位置と色で円を描く †
円を描画する関数
circle=function(x0,y0){
t=seq(0,2*pi,length=100)
x=0.1*cos(t)+x0
y=0.1*sin(t)+y0
polygon(x,y,col=rgb(x0,y0,(x0+y0)/2), lty=0)
}
上の関数をランダムな位置で呼び出す
最終更新：2012年05月28日 20:13
