連立方程式を解く
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
// 蓮見 on Mar. 19, 2012
*パッケージnleqslv [#i28b1d5f]
連立方程式を解くためのパッケージとして、nleqslv(Solve systems of non linear eqations)が用意されている。
同名の関数により、線形連立方程式および非線形連立方程式を解くことができる。
Broyden's methodとニュートン法が選択でき、ヤコビアンをイテレーション毎に計算するかどうかの違いがある。ヤコビアンを自分で与えることもできる。
**関数nleqslv [#dc1d4c53]
例えば、(x,y)を未知数とする二変数の連立方程式を解きたければ、
[f(x, y) = 0, g(x, y) = 0]と変形して、
c(f(x, y), g(x, y))を返す関数を用意する。
関数nleqslvは、c(f(x, y), g(x, y))をゼロにするx, yを求める。
***例: [#w9cd7cd5]
Pareto分布(a, c)の平均と分散は M = a*c/(a-1), V = a*c^2/(a-1)^2/(a-2)だが、
任意の平均と分散を持つPareto分布のパラメータ(a, c)を求める。
fn <- function(x, M, V) {
a <- x[1]
c <- x[2]
c(M - a*c/(a-1), V - a*c^2/(a-1)^2/(a-2))
}
M = 1, V = 1とすると
library(nleqslv)
fn1 <- function(x) fn(x, 1, 1)
ans <- nleqslv(c(3, 3), fn1) # a = 3, c = 3 を初期値とする
結果は
> ans
$x
[1] 2.4142136 0.5857864
$fvec
[1] 1.234427e-10 7.169013e-09
$termcd
[1] 1
$message
[1] "Function criterion near zero"
$scalex
[1] 1 1
$nfcnt
[1] 13
$njcnt
[1] 1
a = 2.4142136, c = 0.5857864が解であることがわかる。
終了行:
// 蓮見 on Mar. 19, 2012
*パッケージnleqslv [#i28b1d5f]
連立方程式を解くためのパッケージとして、nleqslv(Solve systems of non linear eqations)が用意されている。
同名の関数により、線形連立方程式および非線形連立方程式を解くことができる。
Broyden's methodとニュートン法が選択でき、ヤコビアンをイテレーション毎に計算するかどうかの違いがある。ヤコビアンを自分で与えることもできる。
**関数nleqslv [#dc1d4c53]
例えば、(x,y)を未知数とする二変数の連立方程式を解きたければ、
[f(x, y) = 0, g(x, y) = 0]と変形して、
c(f(x, y), g(x, y))を返す関数を用意する。
関数nleqslvは、c(f(x, y), g(x, y))をゼロにするx, yを求める。
***例: [#w9cd7cd5]
Pareto分布(a, c)の平均と分散は M = a*c/(a-1), V = a*c^2/(a-1)^2/(a-2)だが、
任意の平均と分散を持つPareto分布のパラメータ(a, c)を求める。
fn <- function(x, M, V) {
a <- x[1]
c <- x[2]
c(M - a*c/(a-1), V - a*c^2/(a-1)^2/(a-2))
}
M = 1, V = 1とすると
library(nleqslv)
fn1 <- function(x) fn(x, 1, 1)
ans <- nleqslv(c(3, 3), fn1) # a = 3, c = 3 を初期値とする
結果は
> ans
$x
[1] 2.4142136 0.5857864
$fvec
[1] 1.234427e-10 7.169013e-09
$termcd
[1] 1
$message
[1] "Function criterion near zero"
$scalex
[1] 1 1
$nfcnt
[1] 13
$njcnt
[1] 1
a = 2.4142136, c = 0.5857864が解であることがわかる。
ページ名: