読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R {deSolve} で差分方程式

R の deSolve パッケージで差分方程式を計算するには ode のメソッド "iteration" を使う。関数 func は変化率ではなく状態変数の新しい値を返すように書く。 library(deSolve) Ti <- 100 a <- 1.2 disc_logis <-function(Time,x,a){ x2 =a*x*(1-x) list(x2)…

R によるシミュレーテッド・アニーリング法の練習

アルゴリズム 関数 を最大にする を探す問題を最適化問題という. シミュレーテッド・アニーリング法は以下のようなアルゴリズムで最適化問題を解く. 反復 t で, をシミュレートする. 確率 で を受理する. さもなくば とする. ここで . また密度 g は原点 0 …

一次元の拡散方程式をノイマン境界条件の下で解く(差分法、陰解法)

陰解法 以下の偏微分方程式が拡散方程式である。ここで x は空間方向の位置を表し、t は時間を表す変数である。 が求めたい未知関数で、初期条件 は与えられた関数である。D は正の定数である。条件 、 のように端点における微分の値を固定する条…

メモ:第二種不完全ガンマ関数の積分

Wolfram|Alpha: Computational Knowledge Engine #Rで確かめ算 > Igamma <- function(a,x){ + pgamma(x, a, lower=FALSE) * gamma(a) + } > integrate(function(x)Igamma(2,x),3,Inf) 0.2489353 with absolute error < 9.3e-05 > IIgamma <-function(y,k){ +…

ワイブル分布に従う確率変数の和の分布の最尤推定

※初公開時は R のコードが間違っていて、うまくパラメータが求まっていなかった。(修正:7/28) # R code shp1 <- 5 shp2 <- 5 scl1 <- 2 scl2 <- 8 n <-1000 Z <- rweibull(n,shp1,scl1) + rweibull(n,shp2,scl2) LL <- function(par){ m1 <- par[1] eta1 …

R による最急降下法の練習

モチベーション 最尤法とか最小二乗法とかのとき、いまぼくは R の optim 関数をブラックボックス的に使っていて変なとこに収束しちゃうことがよくある。なので最適化について基本的なことを勉強しなおしたい。 追記 「例題」のところにあやまった記述があっ…

メモ:可逆反応の微分方程式

可逆反応の微分方程式 反応物 A が反応物 B を生成し, 生成物から反応物に戻る反応もあるとする.微分方程式,は解,を持つ. R で計算 実線が解析解, 点が数値解. rxnrate=function(t,c,parms){ # rate constant passed through a list called parms k1=parms$k…

微分方程式を含むモデルのパラメータ推定(フェルフルストの人口モデル)

要旨 R で微分方程式を含むモデルのパラメータ推定を行う. まずは単純なモデル(解析的に解ける常微分方程式)で練習. deSolve パッケージの ode 関数(4次ルンゲ=クッタ法)で微分方程式を数値積分しつつ, minpack.lm パッケージの nls.lm 関数(Levenberg…

deSolve パッケージで拡散方程式(熱伝導方程式)を数値的に解く

溶質が均一に分布している材料 A から, 材料 B に溶質が供給されていく状況を考える.一次元の拡散方程式は, ここで, D は拡散係数と呼ばれる定数 (D > 0) である.初期条件と境界条件を以下のように与える.初期条件: 境界条件: この微分方程式を deSolve パ…

deSolve パッケージでロトカ・ヴォルテラの方程式を数値的に解く

捕食者と被食者の増減関係 生物の種は互いに影響しあう.ここでは食う‐食われる(捕食者と被食者;predator-prey)の関係を考える. 被食者は捕食者がいなければ, その個体数に比例して増加する( ) 捕食者は被食者がいなければ, その個体数に比例して減少す…

拡散方程式, 差分法, 陽解法, Rによる練習.

目的 拡散方程式の差分解法として最も簡単と思われる陽解法をとりあえず動かしてみる.※本文中のコード、あってるか自信ない. 拡散方程式の差分表式 一次元の拡散方程式は, ここで, u: 濃度 x: 位置 t: 時間 a: 拡散係数と呼ばれる定数 (a > 0) である.今回は…

オイラー法, テイラー級数法, ルンゲ=クッタ法(R による練習)

オイラー法(Euler's Method) オイラー法(Euler's Method)とは, 1階常微分方程式の数値解法の中でおそらくもっともかんたんなもの. この方法は、数学的に理解しやすく、プログラム的にも簡単なので、数値解析の初歩的な学習問題としてよく取りあげられる…

モンテカルロ法の収束のオーダーについて少し調べた

動機 シミュレーションとか, とりあえず10000回くらい回しとけばいいのかな, くらいの認識しかないので, ちゃんと知っときたい. 定義;オーダー(ランダウの記号) ふたつの関数 と があるとする. と表す(これを f が のオーダーであるという)とき,その意…

感染症のモデル(SIRモデル)に入門した

導入 記法はウィキペディアに合わせる. SIRモデル - Wikipedia 時刻 t において, S(t):感染可能者(Susceptible)の数. これから病気にかかるおそれのある人たち. I(t):感染者(Infected)の数. いま病気にかかっていて人に病気を移す可能性がある人たち. …

一次元の解探索、f(x)=0を解く、二分法(bisection method)

大雑把にいうと、 f(a)とf(b)の符号が違えば、aとbの間にf(c)=0となるようなcがあるはずだ だから区間[a,b]を半分づつ狭めていってcに近づこう というのが二分法の考え方である。 二分法(wikipedia) なんてわかりやすいアルゴリズムなんだ こういうのがい…

一次元の解探索、f(x)=0を解く、uniroot()、ブレント法(Brent's method)

(前回) いま読んでいるMaria L. Rizzo "Statistical Computing with R"では, 「ブレント法(Brent's method)は,root bracketingと二分法(bisection),逆二次補間の組み合わせである」 と書いてあるんですが,ウィキペディアでは 「ブレント法は、二分…

2進法、コンピュータ、小数、恒等、ニアリーイコール

Statistical Computing with R (Chapman & Hall/CRC The R Series)作者: Maria L. Rizzo出版社/メーカー: Chapman and Hall/CRC発売日: 2007/11/15メディア: ハードカバー クリック: 14回この商品を含むブログ (6件) を見るちょっと煮詰まってきたので仕方な…

オーバーフロー、アンダーフロー

Statistical Computing with R (Chapman & Hall/CRC The R Series)作者: Maria L. Rizzo出版社/メーカー: Chapman and Hall/CRC発売日: 2007/11/15メディア: ハードカバー クリック: 14回この商品を含むブログ (6件) を見る(前回) Maria L. Rizzo "Statist…

integrate()で数値積分

統計的データ解析のための数値計算法入門 (統計ライブラリー)作者: 岩崎学出版社/メーカー: 朝倉書店発売日: 2004/10メディア: 単行本購入: 1人 クリック: 2回この商品を含むブログ (5件) を見るこの本の例5.1.1より、 確率変数Xの確率密度関数が のときの、…