2010年12月21日

科研費の繰越が可能に

科学研究費、繰り越し可能に…一部を基金で運用
http://headlines.yahoo.co.jp/hl?a=20101221-00000172-yom-sci
読売新聞 12月21日(火)8時6分配信

引用開始) 政府は、来年度予算編成で、「科学研究費補助金(科研費)」の一部について、翌年度への繰り越しが自由にできる制度改正を行う方針を決めた。

 単年度予算の制約を受けない「基金」の形で補助金を運用する。年度末の予算消化のために生まれる無駄をなくし、研究費を効率的に使える体制に改める。研究が進んだ場合は、研究費を前倒しして使うこともできる。

 1965年度に創設された科研費は現在、年間約6万件の研究を支える日本の学術研究の土台だが、抜本的な制度改正は初めて。

 総額2000億円(今年度)に上る科研費のうち、来年度はまず約310億円を基金化する。若手研究者向けの研究費が中心になる。現在でも研究計画に変更があった場合は国に申請すれば繰り越しが認められるが、そのための時間や労力が研究の妨げになるとして、研究者が長年、制度改正を強く求めていた。(引用終わり

いやはや、僕はあまりお金使わなくても出来る研究なので科研費と結構無縁に学生生活送ってきてたのですが、横国に移ってからは科研費を使ってる人がとっても多くてこのニュースはありがたいんではないかなぁ、と思いました。実際、知り合いもあと○万円使わないといけないんですよー、とか言ってたし、使いきらなくても減らされないならその方が良いのだろう。

posted by しばきん at 16:15| 神奈川 ☁| 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) | 数学 | このブログの読者になる | 更新情報をチェックする

逆行列の勉強

行列式、正直言って苦手なのです。何やってるかイメージわかない。でも必要になってきたので、お勉強します。とりあえず逆行列ってどんなときに使うのか勉強しました。

参考元URL

> #######################
> ###逆行列はどんなときに使うのか###
> #######################
>
> ###行列の基本性質のおさらい###
> ###ベクトルを合わせたものとしての行列###
> A <- as.matrix(cbind(c(10, 2, 3), c(4, 50, 6), c(7, 8, 90)))
> b1 <- as.matrix(c(100, 11, 12))
> b2 <- as.matrix(c(13, 140, 15))
> b3 <- as.matrix(c(16, 17, 180))
>
> ###これで行列の積となってる###
> A%*%b1
[,1]
[1,] 1128
[2,] 846
[3,] 1446
>
> ###ベクトルだと思ってチコチコ"縦のまま"計算###
> ###今更ながらに重回帰のときの行列表記にしっくり感###
> A[, 1]*b1[1, 1]+A[, 2]*b1[2, 1]+A[, 3]*b1[3, 1]
[1] 1128 846 1446
>
> ###他のも計算###
> A[, 1]*b2[1, 1]+A[, 2]*b2[2, 1]+A[, 3]*b2[3, 1]
[1] 795 7146 2229
> A[, 1]*b3[1, 1]+A[, 2]*b3[2, 1]+A[, 3]*b3[3, 1]
[1] 1488 2322 16350
>
> ###普通に行列計算したって一緒###
> ###なるほどAは生データbは係数だと考えるとすごく納得###
> B <- cbind(b1, b2, b3)
>
> ###係数を取りだすという意味からの逆行列###
> ###逆行列はその性質から元の行列と掛け合わせると単位行列となる###
> ###なんか整数がいいので意味一緒だけどroundで###
> round(solve(A)%*%A)
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
>
> ###これでも同じ###
> round(cbind(solve(A)%*%A[, 1], solve(A)%*%A[, 2], solve(A)%*%A[, 3]))
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
>
> ###さらにA列にそれぞれ何か値が掛かっているとする###
> A2 <- A %*% c(1, -2, 30)
>
> ###先ほどの逆行列を掛けてやると掛けていた値が得られる###
> solve(A)%*%A2
[,1]
[1,] 1
[2,] -2
[3,] 30
>
> ###何故かというと線形性があるから###
> ###つまり線形和として下の式のかけ算の後ろの項を見たときにその係数を取りだす作業になっている###
> solve(A)%*%(A[, 1]*1 + A[, 2]*-2 + A[, 3]*30)
[,1]
[1,] 1
[2,] -2
[3,] 30
>
> ###行列がかかっているなら###
> A3 <- A%*%as.matrix(cbind(c(1, -2, 30), c(40, 5, -6), c(-7, 80, 9)))
> solve(A)%*%A3
[,1] [,2] [,3]
[1,] 1 40 -7
[2,] -2 5 80
[3,] 30 -6 9
>
> ###対角化のための逆行列###
> ###行列の固有値と固有ベクトル###
> eigen(A)
$values
[1] 91.471925 48.937866 9.590208

$vectors
[,1] [,2] [,3]
[1,] 0.09339418 -0.07440825 -0.99850799
[2,] 0.19292236 -0.98595733 0.04267442
[3,] 0.97675918 0.14950438 0.03406895

>
> ###固有値のみを対角成分で抽出###
> round(solve(eigen(A)$vectors)%*%A%*%eigen(A)$vectors, 2)
[,1] [,2] [,3]
[1,] 91.47 0.00 0.00
[2,] 0.00 48.94 0.00
[3,] 0.00 0.00 9.59
> #######################
> ###逆行列はどんなときに使うのか###
> #######################
ラベル:R 数学
posted by しばきん at 01:13| 神奈川 ☀| Comment(0) | TrackBack(0) | 数学 | このブログの読者になる | 更新情報をチェックする

2010年12月15日

自作通常クリギングプログラム

行列でパラメータ推定することを覚えたので、クリギングを自作プログラミングしてみました。なるほど、いつもクリギングは時間がかかるんだなぁと思ってましたが、こういう計算をしているのであれば納得です。

参考文献
クリギングとその地理的応用

> ####################
> ###自作クリギングプログラム###
> ####################
> #pdf("kriging.pdf")
> set.seed(2)
> library(lattice)
> library(spatstat)
> library(gstat)
> ###生データの作成###
> raw <- data.frame(x=rnorm(100, 0, 1), y=rnorm(100, 0, 1), z=0)
> raw[, 3] <- log(10*abs(1/(raw[, 1]*raw[, 2])))
> raw2 <- raw
>
> ###生データのサンプル場所の確認###
> #plot(raw2[, 1], raw2[, 2])
>
>
> ###まずは既存のパッケージを使ってやってみる###
>
> ###座標系としてやる###
> coordinates(raw2) = ~ x + y
>
> ###適当にあてはめ###
> v <- variogram(z~1, raw2)
> #plot(v)
>
> ###初期値を与えて当てはめ###
> v.fit <- fit.variogram(v, vgm(psill=0.9, "Gau", range=7, nugget=0.2))
>
> ###シルとかレンジの値###
> #summary(v.fit)
> v.fit
model psill range
1 Nug 0.1704907 0.0000000
2 Gau 1.8520082 0.6434817
>
> ###点に含まれる数が等しくないと当てはまりが悪いように見えることもある###
> #plot(v, model=v.fit)
>
> ###予測グリッドの作成###
> mesh.grid <- expand.grid(x=seq(-3, 3, 0.1), y=seq(-3, 3, 0.1))
>
> ###クリギング補間###
> gs <- gstat(id="z", formula=z~1, data=raw2, model=v.fit)
> coordinates(mesh.grid) = ~x+y
> res_gs <- predict.gstat(gs,mesh.grid)
[using ordinary kriging]
> aa <- data.frame(z=res_gs$z.pred ,expand.grid(x=seq(-3, 3, 0.1),
+ y=seq(-3, 3, 0.1)))
> ###推定された図###
> levelplot(aa[, 1] ~ aa[, 2] + aa[, 3], main="package_used")
>
> ###平均値と分散を確認###
> mean(aa[, 1])
[1] 2.661233
> var(aa[, 1])
[1] 0.7555209
>
>
> ###次に自作でプログラムを組んで結果が一緒になるか確認###
>
> ###地点間の距離と分散を計算###
> d_temp <- u_temp <- v_temp <- a <- a2 <- NULL
> for(i in 1 : nrow(raw)) {
+ d_temp <- sqrt((raw[i, 1]-raw[, 1])^2+(raw[i, 2]-raw[, 2])^2)
+ u_temp <- (raw[i, 3]+raw[, 3])/2
+ v_temp <- (raw[i, 3]-u_temp)^2+(raw[, 3]-u_temp)^2
+ a <- cbind(d_temp, v_temp)
+ a2 <- rbind(a2, a)
+ }
>
> ###バリオグラムにモデルを当てはめるときのカットポイントを指定###
> ###ここが違うので既存パッケージの結果とぴったり一緒のバリオグラムモデルは推定されないことに注意###
> bp <- seq(0, 2.1, length=16)
>
> ###切った間の平均値を計算してそれにモデルを当てはめる###
> a2_temp <- a3 <- a4 <- NULL
> for(i in 1 : (length(bp)-1)){
+ a2_temp <- subset(a2, a2[, 1]<=bp[i+1]&a2[, 1]>=bp[i])
+ a3 <- apply(a2_temp, 2, mean)
+ a4 <- rbind(a4, a3)
+ }
>
> ###確認###
> #plot(a4, pch=16)
>
> ###ガウスモデルのパラメータ推定###
> func_vario <- function(par, data, dis) {
+ inter <- par[1]
+ phi1 <- par[2]
+ phi2 <- par[3]
+ return(
+ sum(
+ (data-(inter + phi1*(1-exp(-(dis/phi2)^2))))^2
+ )
+ )
+ }
> res <- optim(c(1, 2, 1), func_vario, data=a4[, 2], dis=a4[, 1], hessian=T)
>
> ###出力用に関数を作る###
> func_vario_res <- function(inter, phi1, phi2, dis) {
+ return(
+ inter + phi1*(1-exp(-(dis/phi2)^2))
+ )
+ }
>
> ###さっきのプロットに線を重ねる###
> #points(a4[, 1], func_vario_res(inter=res$par[1], phi1=res$par[2],
> # phi2=res$par[3], dis=a4[, 1]), type="l", lwd=2
> #)
>
> ###任意の地点間の距離を計算してくれる関数の作成###
> func_dist <- function(pred, obs) {
+ pred_p <- ppp(pred[, 1], pred[, 2], c(-1e+5, 1e+5), c(-1e+5, 1e+5))
+ obs_p <- ppp(obs[, 1], obs[, 2], c(-1e+5, 1e+5), c(-1e+5, 1e+5))
+ return(crossdist(obs_p, pred_p))
+ }
>
> ###生データ間の距離を計算###
> dis <- func_dist(pred=raw, obs=raw)
>
> ###距離を使ってバリオグラムの計算###
> X <- func_vario_res(inter=res$par[1], phi1=res$par[2],
+ phi2=res$par[3], dis=dis)
>
> ###重みの合計は1###
> X <- rbind(X, 1)
> X <- cbind(X, 1)
>
> ###対角成分は0距離###
> diag(X) <- 0
>
> ###予測グリッドの作成###
> mesh.grid <- expand.grid(x=seq(-3, 3, 0.1), y=seq(-3, 3, 0.1))
>
> ###予測グリッドと生データの距離を計算###
> dis2 <- func_dist(pred=mesh.grid, obs=raw)
>
> ###バリオグラムの計算###
> for(i in 1 : nrow(mesh.grid)){
+ X2 <- func_vario_res(inter=res$par[1], phi1=res$par[2],
+ phi2=res$par[3], dis=dis2[, i])
+ X2 <- c(X2, 1)
+
+ ###クリギング方程式を解く###
+ ###一番下はラグランジュの未定係数法のパラメータであって推定には使わない###
+ mesh.grid[i, 3] <-sum(raw[, 3]*as.numeric(X2%*%solve(X))[1:100])
+ }
> #x11()
> levelplot(mesh.grid[, 3] ~ mesh.grid[, 1] + mesh.grid[, 2], main="original_program")
>
> ###パッケージで推定したものとほとんど同じ平均値と分散###
> mean(mesh.grid[, 3])
[1] 2.662301
> var(mesh.grid[, 3])
[1] 0.7765148

> ###4分位点も確認。大体一緒###
> quantile(aa[, 1])
0% 25% 50% 75% 100%
1.028423 2.156357 2.553124 2.760064 6.581620
> quantile(mesh.grid[, 3])
0% 25% 50% 75% 100%
0.974662 2.154865 2.559820 2.748471 6.761476
> #dev.off()
> ####################
> ###自作クリギングプログラム###
> ####################


出来た図。kriging.pdf
パッケージを使ったものとほとんど同じ結果。なぜ違うかというと、バリオグラムのパラメータ推定時のカットポイントの値が違うから(だと思う)。

posted by しばきん at 20:23| 神奈川 ☁| Comment(0) | TrackBack(0) | 空間統計 | このブログの読者になる | 更新情報をチェックする

行列式による回帰係数の導出

そういえば、最小二乗法で回帰分析するときにはラクしてlm関数ばっかり使ってるけど、行列式でもちゃんと理解しておこう、と思って随分前にダウンロードした以下の資料でお勉強。
行列の転置と微分

ちゃんと一致するか確認。


> set.seed(1)
> a <- 2
> b <- 3
> x <-seq(1, 10)
> y <- a + b*x + rnorm(10)
###切片の分を入れてやる###
> x2 <- matrix(c(rep(1, 10), x), 10, 2)
> solve(t(x2)%*%x2)%*%t(x2)%*%y 
[,1]
[1,] 1.831176
[2,] 3.054732
> coef(lm(y~x))
(Intercept) x
1.831176 3.054732 #当然一致する。




ラベル:R 数学
posted by しばきん at 16:50| 神奈川 ☁| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする
×

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