2012年03月19日

黄金比率とゼロと無限

今日読んだ本に触発されて、唐突にフィボナッチ数列で遊んでみたくなった。
いつものごとくRで。

> Fibo <- function (n) { 
+ phi <- (1 + sqrt(5))/2
+ return(
+ (phi^n - (-phi)^-n)/sqrt(5)
+ )
+ }
>
> Fibo(1:10) #実験
[1] 1 1 2 3 5 8 13 21 34 55
>
> x <- seq(1, 30)
> Fibo(x+2)/Fibo(x+1) #黄金比に収束
[1] 2.000000 1.500000 1.666667 1.600000 1.625000 1.615385 1.619048 1.617647
[9] 1.618182 1.617978 1.618056 1.618026 1.618037 1.618033 1.618034 1.618034
[17] 1.618034 1.618034 1.618034 1.618034 1.618034 1.618034 1.618034 1.618034
[25] 1.618034 1.618034 1.618034 1.618034 1.618034 1.618034

おおー。ついでに対数螺旋に勝手にz軸をこさえて遊んでみる。

> library(aspace)
> library(scatterplot3d)
>
> func <- function(theta, a, b){
+ a*exp(b*theta)
+ }
> degree <- seq(-1e+4, 1e+2)
>
> x <- func(theta=degree/180*pi, a=0.1, b=0.2) * cos(degree/180*pi)
> y <- func(theta=degree/180*pi, a=0.1, b=0.2) * sin(degree/180*pi)
> jpeg("Golden_number.jpeg")
> plot(x, y, type="l", col=4, lwd=3)
> dev.off()
> jpeg("Golden_number_0_8.jpeg")
> scatterplot3d(1/(x*y) ~ x + y, type="l", color=4, zlim=c(-1e+30, 1e+30))
> dev.off()

対数螺旋
Golden_number.jpeg

対数螺旋とz軸
Golden_number_0_8.jpeg
うーん、マイナス∞と∞が座標(0, 0)を通る垂直線上に現れてくる…。あぁすごいなぁ。
ラベル:数学 R
posted by しばきん at 23:10| 神奈川 ☁| Comment(0) | TrackBack(0) | 数学 | このブログの読者になる | 更新情報をチェックする

2011年06月16日

動的計画法による最適漁獲方策の導出練習

昨日の線形計画法に続き次は動的計画法を勉強しました。
逆瀬川浩孝さんのプログラムを見ながら勉強しました。ありがとうございます。
この順番で学習するのが最適なのかは分かりません。以下有名な(?)ナップザック問題を水産風にアレンジしたものです。・・・といっても別に名称を変えただけです。練習練習。

###漁師が漁に行く場合を考える###
###漁獲物を保管する場所は一つしかなく、容量も決まっているとする###
###重い大きな魚ほど価値が高いとする###
###従って容量が決まっている時に、どれほどの割合で大きな魚とそうでない魚を漁獲すれば###
###最も効率が良いのかを解く###
###なお魚は4種が1尾づついるとし、重複して同じ種は漁獲出来ないものとする###

###漁獲物の重量を決める###
> w <- c(4,6,3,2)
>
> ###漁獲物の価値を決める###
> v <- c(13,18,12,8)
>
> ###計算結果を格納する行列を用意する###
> ###1行目と1列目は計算の便宜上使うだけ。最後は捨てる###
> M <- matrix(0, length(w)+1, sum(w)+1)
>
> ###漁獲物を入れるスペース###
> T <- sum(w)+1
>
> ###ある一つの漁獲物ごとにスペースに入るのか否かまず計算する###
> for (n in 1 : length(w)) {
+ for (t in 1 : (T-1)) {
+ if(w[n] > t) {
+
+ ###入らなければ一つ上の値を持ってくる###
+ ###2つめの漁獲物が入らない時は、既に入っている1つめの漁獲物の値が挿入される###
+ M[n+1, t+1] <- M[n, t+1]
+ } else {
+
+ ###もし入るなら1つ上の値と、2つ目の漁獲物の価値と2つ目が入ってなお残るスペースで入っている漁獲物###
+ ###の合計を比べて大きい方を採用###
+ ###こうすることで別に軽い順に入れていかなくてもあとで更新される仕組みとなっている###
+ M[n+1, t+1] <- max(M[n, t+1], v[n] + M[n, t-w[n]+1])
+ }
+ }
+
+ }
>
> ###余計な計算行列を削除して最終出力###
> ###行の名前に重さと価値を小数点で区切って両方載せてみる###
> M2 <- M[(2:nrow(M)), (2:ncol(M))]
> rownames(M2) <- v+0.1*w
> rownames(M) <- c(0, v+0.1*w)
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
13.4 0 0 0 0 13 13 13 13 13 13 13 13 13 13 13 13
18.6 0 0 0 0 13 13 18 18 18 18 31 31 31 31 31 31
12.3 0 0 0 12 13 13 18 25 25 30 31 31 31 43 43 43
8.2 0 0 8 12 13 20 21 25 26 33 33 38 39 43 43 51
> M2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
13.4 0 0 0 13 13 13 13 13 13 13 13 13 13 13 13
18.6 0 0 0 13 13 18 18 18 18 31 31 31 31 31 31
12.3 0 0 12 13 13 18 25 25 30 31 31 31 43 43 43
8.2 0 8 12 13 20 21 25 26 33 33 38 39 43 43 51

上記Mがどういうメカニズムで埋められたのか日本語に起こしてみる。ここでの注意は、計算上作った1行目と1列目のせいで、例えばn=2であるとしたとき、重さ-価値関数(w[2], v[2])はちゃんと2番目の漁獲物の重さと価値を引いてくるが、計算結果を入れる行列Mの中ではM[2, ]としたときは漁獲物1を意味しているということである。以下順番に追ってみる。
例えばMの7列目(t=6)を見ると、
    [,7]
0 0
13.4 13
18.6 18
12.3 18
8.2 21

となっている。まずn=1のときに、最初に重さ4、価値13の魚が2行目に入った。次にn=2のとき、つまり3行目に重さ6の魚が入る際、1つ上の重さ4の魚の価値と、(入るならば)1つ上の重さ4の魚と新しく重さ6の魚を足し合わせた価値を比べる。つまり、今容量が6(t=6)あるわけだから6-4+1(0列目分)=3(t=2のときの状態)が余っている容量であるので3列目に戻る。
    [,6-4+1]
0 0
13.4 0
18.6 0
12.3 0
8.2 8

ここでは、2行目(n=2)に値が入っていない。容量が6-4=2のときには、重さ4の魚は入れないからである(プログラム中では+1されて一つ横の列に挿入される)。従って、重さ6の魚と重さ4の魚は、容量が6しかないときには同時には入れることが出来ない、と判定される。だからM2の7列目の3行目(M[2+1,6+1])には(つまりn=2, t=6)、1つ上の値(M[2,6+1]、つまりn=2,t=6)13と、18(v[2])+0(M[2,6-4+1])=18という数字を比べて、大きい18という数字が挿入される。どちらか1つしか入らないのであれば、重さ4の魚と重さ6の魚を比べて価値の大きい方を取った、というわけである。n=3のとき、つまり4行目に値を入れる場合も同様。容量6-3+1=4列目は、
     [,6-3+1]
0 0
13.4 0
18.6 0
12.3 12
8.2 12

なのでM[3+1, 6+1]には、前の状態までで最適であったM[3, 6+1]=18と、12(v[3])+0(M[3,6-3+1])=12とを比べて、今回は更新されず18のままとなる(M[4,7])。最後に、n=4つまり5行目に値を入れる。容量6-2+1列目は、
     [,6-2+1]
0 0
13.4 13
18.6 13
12.3 13
8.2 13

なのでM[4+1, 6+1]には、前の状態までで最適であったM[4, 6+1]=18と、8(v[4])+13(M[4,6-2+1])=21を比べて、今度は21が採択されるのである。こうして、容量が6のときには、価値が8のものと13のものを組み合わせるのがベストだ、ということが分かりました。重さを引いていけば、どういう組み合わせで入ったのかは逆算できます。前の状態までで最適であったものを使うのがこの動的計画法のミソである。・・・と理解しました。

まだ自分一人で使えるレベルにないので、引き続き練習します。
ラベル:R 水産資源 数学
posted by しばきん at 17:16| 神奈川 ☁| Comment(0) | TrackBack(0) | 数学 | このブログの読者になる | 更新情報をチェックする

2011年06月15日

漁業者の最適漁獲方策を線形計画法で解いてみる練習

数理計画法の勉強のためにまずは線形計画法を勉強しました。それが最適な勉強法なのかはわかりません。水産出身なので、漁業者を例にとります。あと、Seesaで埋め込み数式をlatex形式で書けるらしいので、その勉強の意も込めて練習練習。

・漁業者が一人いて2魚種を漁獲する場合を考える
・魚種ごとにかかるコストも利益も違うとした
・どちらの魚種をそれぞれどれほど漁獲すれば最も利益が上がるか解く
・コストには漁場に往復にかかるコストと着いてから操業にかかるコストの2つがあるとした

まず魚種は2魚種なので、それぞれの資源量を、コストは往復にかかるコストと操業にかかる合計のコストをそれぞれ、として、それらがそれぞれの魚種にだけかかるものとする。利益としてそれぞれの魚種をだけ漁獲した際の利益をとする。

・最大化する項


・制限する項(資源量の10%は置いておく)





・非負条件




では、これをRに解かせてみます。
> ###2魚種の個体数###
> X1 <- 10
> X2 <- 100
>
> ###利益###
> p1 <- 1
> p2 <- 0.5
>
> ###コスト(往復)###
> c1 <- 1
> c2 <- 0.8
> C <- 100
>
> ###コスト(操業)###
> c3 <- 1
> c4 <- 0.2
> C2 <- 50
>
> ###最大化する項はp1*x1+p2*x2###
> ###pの値だけ渡す###
> Profit <- c(p1, p2)
>
> ###制約の上限###
> ###最近WormとHilbornの研究を読んだのにちなんで###
> ###バイオマスの10%を下回らないようにする###
> Limit <- c(C, C2, X1*0.9, X2*0.9)
>
> ###制限する項###
> ###それぞれの個体数も制限条件に入る###
> Cost1 <- c(c1, c2)
> Cost2 <- c(c3, c4)
> N1 <- c(1, 0)
> N2 <- c(0, 1)
> Cost <- rbind(Cost1, Cost2, N1, N2)
>
> ###解###
> solveLP(Profit, Limit, Cost, TRUE)
> S <- solveLP(Profit, Limit, Cost, TRUE)$solution
>
> ###X1を十分大きい数にして再試行###
> X1 <- 100
> Limit <- c(C, C2, X1*0.9, X2*0.9)
> S2 <- solveLP(Profit, Limit, Cost, TRUE)$solution
>
> ###図示して確認###
> pdf("LP.pdf")
> plot(0, xlim=c(0, 100), ylim=c(0, 100), col=0, ylab="x1", xlab="x2")
> abline(a=C, b=-c2, lwd=2, col=2)
> abline(a=C2, b=-c4, lwd=2, col=4)
> points(S[2], S[1], col=1, cex=2, pch=16)
> points(S2[2], S2[1], col=3, cex=2, pch=16)
>
> legend(60, 100, c("Time_cost", "Operating_cost"), lwd=2, col=c(2, 4))
> legend(60, 85, c("Optimal_catch_X1_High"), pch=16, col=3)
> legend(60, 75, c("Optimal_catch_X1_low"), pch=16, col=1)
> dev.off()


出来た図。
LP.pdf
資源量が少ない時(黒点)は、コストから導かれる最適な漁獲が出来ずそれを下回ったところで最適解となることが分かります。当たり前ですけれども、まぁ練習なので。
ラベル:数学 水産資源 R
posted by しばきん at 12:38| 神奈川 ☁| Comment(0) | TrackBack(0) | 数学 | このブログの読者になる | 更新情報をチェックする

2011年01月29日

数学の記号

大変有難い。
数学の記号まとめ

「属する(ぞくする)」っていうのを単なる和の意味だと思って本を読んでたら最後になって計算が合わなくなってえらい目にあったので、自戒を込めて。
ラベル:数学
posted by しばきん at 00:34| 神奈川 ☀| Comment(0) | TrackBack(0) | 数学 | このブログの読者になる | 更新情報をチェックする

2010年12月18日

逆行列の意味=割り算

昨日に引き続き、逆行列の意味を考えてみた。割り算だと考えることもできそう。
参考文献URL1
参考文献URL2

#######################################
> ###回帰係数を眺めていて気付いた逆行列の意味=割り算###
> #######################################
> set.seed(1)
> a <- 2
> b <- 3
> n <- 100
> x <-seq(1, n)
> y <- a + b*x + rnorm(n)
> ##切片の分を入れてやる###
> x2 <- matrix(c(rep(1, n), x), n, 2)
>
> ###行列で###
> solve(t(x2)%*%x2)%*%t(x2)%*%y 
[,1]
[1,] 2.131666
[2,] 2.999549
>
> ###パッケージで###
> coef(lm(y~x))
(Intercept) x
2.131666 2.999549
>
> ###最小二乗推定量から###
> b2 <- sum((x2[, 2]-mean(x2[, 2]))*(y-mean(y)))/sum((x2[, 2]-mean(x2[, 2]))^2)
> a2 <- mean(y)-b2*mean(x2[, 2])
> a2;b2
[1] 2.131666
[1] 2.999549
>
> ###どうせ分母分子n-1で割るんだからvarでもOK###
> b3 <- var(x2[, 2], y)/var(x2[, 2])
> a3 <- mean(y)-b3*mean(x2[, 2])
> a3;b3
[1] 2.131666
[1] 2.999549

ここで気づいたのはsolve(t(x2)%*%x2)とsum((x2[, 2]-mean(x2[, 2]))^2)
が同じ意味なのではないか、ということです。確かに逆行列をかけるということは
以前にも見たように係数を取りだすとか、対角化を行うとか色々ありましたが、
それってつまり割り算をしているからそういう風になるんじゃ?
では他に割り算をするようなケースってどんな場合?
じゃあ例えば相関係数でも。

> ###まずは普通に###
> var(x, y)/(sd(x)*sd(y))
[1] 0.9999467
>
> ###パッケージで確認###
> cor(x, y)
[1] 0.9999467
>
> ###何も参考文献なしで行列化を試みる###
> solve(sqrt(t(x-mean(x))%*%(x-mean(x))/(n))%*%
+ sqrt(t(y-mean(y))%*%(y-mean(y))/(n)))%*%
+ t(x-mean(x))%*%(y-mean(y))/(n)
[,1]
[1,] 0.9999467

出来た!割り算だと思うと、逆行列の見方も変わりますなぁ。
一応参考文献見てそれっぽいやり方もやってみる。

> ###参考文献見て書いた相関行列の算出法###
> ###基本的にやっていることは一緒ということを確認###
> ###あらかじめ平均値を引いておく###
> Xc <- as.matrix(cbind(x-mean(x), y-mean(y)))
>
> ###分散共分散行列を作成###
> Sxx <- 1/n*t(Xc)%*%Xc
>
> ###対角成分に入っている分散だけ取り出す###
> Ds <- diag(2)
> diag(Ds) <- diag(Sxx)
>
> ###分散にルートをかけて標準偏差に直してそれで平均値を引いた値をそれぞれ割る###
> Z <- Xc%*%solve(sqrt(Ds))
>
> ###その値を乗じれば相関係数となる###
> Rxx <- 1/n*t(Z)%*%Z
> Rxx
[,1] [,2]
[1,] 1.0000000 0.9999467
[2,] 0.9999467 1.0000000
> #######################################
> ###回帰係数を眺めていて気付いた逆行列の意味=割り算###
> #######################################

でもきっと他の解釈もあるはず。なんか「変換」がキーワードなのかな?引き続き勉強しよう。
ラベル:R 数学 統計学
posted by しばきん at 17:57| 神奈川 ☀| Comment(0) | TrackBack(0) | 数学 | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。