ニュートン法により、方程式の実根を求める。
xn = xn-1 - f(xn-1)/f'(xn-1) (f'=df/dx)
を繰り返して根に近づく。参照：「方程式の根」
x0は初期値。その値により異なる解に達したり、解が得られないことがある。
トレース版では、表示場所を次のようにする。
メイン： xFxNewtonTrace("〇〇", fx, x0, eps);
ＨＴＭＬ： <div id="〇〇">ここに表示される</div>
a3を省略すると２次方程式、さらにa2を省略すると１次方程式になる。
２次方程式は根の公式による。３次方程式は、特殊なケースを除き、ニュートン法により一つの実根を求め、式の変形により２次方程式にして求める。
xfx3(2, 1) x + 2 = 0 １次方程式
xfx3(1, 1, 1) x2 + x + = 0 ２次方程式 根の公式
xfx3(0, 1,-1, 1) x(x2 -x1 +1) = 0 ( ）内は２次方程式
xfx3(2, 0, 0, 1) x3 = -2 ３乗根の公式を使う
xfx3(2, 0,-1, 1) x3 - x2 + 2 = 0 一般的な３次方程式
戻り値は２次元配列である。
var ary = xfx3(～);
とすると、
ary[i][0] 解ｉの実数部、ary[i][0] 虚数部
が入る。
数式 f(x) が最小となるｘの値を戻します。
xmin～xmaxを等分してdxとし、for (x=xmin, x<=xmax; x=x+dx) でf(x) を最小にする x=xpot を求め、次に xopt-dx～xopt+dx で同様な操作を繰り返すことにより、最終的な xopt を求めます。epsを省略すると0.0001になります。
数式 f(x) が最大となるｘの値を戻します。
xmin～xmaxを等分してdxとし、for (x=xmin, x<=xmax; x=x+dx) でf(x) を最大にする x=xpot を求め、次に xopt-dx～xopt+dx で同様な操作を繰り返すことにより、最終的な xopt を求めます。epsを省略すると0.0001になります。
fxのｘにおける微係数aを求めます。ｘの近傍の２点 x-d*xとx-d*xでの勾配をa0とし、ｄを1/2にしたときの勾配aとの差が、相対誤差 |(a-a0)/a| < eps になるまで繰り返します。
dを省略すると0.01、epsを省略すると0.001にしています。
シンプソンの公式（参照：「定積分」）により計算しています。
最初のn2で計算した値をx0、次にn2を２倍（すなわちキザミ幅を1/2）にして計算した値をｘとして、相対誤差 |(x-x0)/x|
n+1組の点xa[0]～xa[n]、ya[0]～ya[n] を通るｎ次式をラグランジュの補間法により求め、m 個の、xa[n+1]～xa[n+1} における ya[n+1]～ya[n+m] を戻します。
４(=n+1)点(0,1), (2,5), (3,8), (5,10) を通る３(=n)次式により、
２(=m)点、ya[4] = y(4), ya[5] = y(6) を求める
var n = 3;
var m = 2;
var xa = [0, 2, 3, 5, 4, 6];
var ya = [1, 5, 8, 10];
polArray(n, m, xa, ya);
解： ya[4] = y(4) = 10.1 ya[5] = y(6) = 6.6
n+1個の点、(x0,y0),(x1,y1),(x2,y2)～(xn,yn) を通るｎ次式をラグランジュの補間法により求め、ｘにおけるｎ次式の値を戻す。
y = polPoints(2, 1,2, 3,6) (1,2)と(3,6) を通る１次式の x=2 の値が y に入る
y = polPoints(2, 0,0, 2,0, 3,3) ３点 (0,0).(2,0),(3,3)を通る２次式の x=-1 値が y に入る
p組の(x,y)から、最小二乗法により、ｎ次の近似式 y=a0+a1x+a2x2+...+anxn の係数 ai を求め。この近似式により、ｆ個の x に対応する y を計算します。x,y,a は、それぞれ一次元配列 xa, ya, aa でメイン側で定義します。
参照：「回帰分析・最小二乗法」、
「多項式、高次式での最小二乗法」
配列 a[i][j] を、j=keycol の列の大小によりソートします。
desc：省略時あるいは正の値なら昇順、負の値なら降順
fixrow：省略時あるいは０ならすべての行を対象。先頭のｎ行を固定する（ソートの対象にしない）ときはｎの値
（fixrowを与えるときは、descは省略できない）
入力データ
a[0] = ["C", 30];
a[1] = ["B", 20];
a[2] = ["E", 50];
a[3] = ["A", 10];
a[4] = ["D", 40];
keycol = 1;
配列ａについて、指定列keycolの値により昇順または降順にソートしたときの各行の順位を戻します。行列そのものはソートされません。
順位[k] = i とは、昇順ソートしたとき、k番目になる行が i であることを示します。
すなわち、行列[[順位[k]][j] がｋ番目の行になります。
a：配列
氏名 学部 学年
a[0][0] = '阿部', a[0][1] = "工学部", a[0][2] = 3
a[1][0] = '井上', a[1][1] = "文学部", a[1][2] = 2
a[2][0] = '宇田', a[2][1] = "文学部", a[2][2] = 1
a[3][0] = '江川', a[3][1] = "工学部", a[3][2] = 4
a[4][0] = '大西', a[4][1] = "法学部", a[4][2] = 3
keycol：順位を調べる列。学年ならば２
desc：省略あるいは正数ならば昇順、負数ならば降順
fixrow：先頭の固定（比較の対象にしない）行数
省略すると全行が比較対象になる。２とすると、先頭の２行（阿部と井上）は比較対象外
ary = orderMatrix(a, 2); とすると、学年で昇順にソートされるので、
ary(0)：０番目に小さい行→宇田（２）→ary(0) = 2
ary(1)：１番目に小さい行→井上（１）→ary(1) = 1
：
ary(4)：４番目に小さい行→江川（３）→ary(3) = 3
となります（同じ値の場合は添え字の小さいほうが先になります）。
aryのサイズは、０～行数－fixrowになります。
昇順に取り出したいのであれば、
for (k=0; k
alert( a[i][0] );
}
とすればよいのです。
入力データ
ａはサイズｎ×ｎの行列、ｂはサイズｎの定数項で、ａｘ＝ｂの解ｘを配列aryに戻します。ａ．ｂのサイズが大きいときは先頭のｎ個が計算対象になります。
ｎを省略するとｂのサイズになります。
入力配列 a と b は、関数呼出後も値が変わりません。
アルゴリズムは「連立方程式」を参照。誤差の考慮はしていないので、ｎが大きな問題は解けません。
function GaussJordanTest1() {
var n = 3;
var a = new Array();
a[0] = [0, 3, 2]; // 0x0 + 3x1 + 2x2 = 12
a[1] = [2, 4, 0]; // 2x0 + 4x1 + 0x2 = 10
a[2] = [3, 1, 5]; // 3x0 + 1x1 + 5x2 = 20
var b = [12, 10, 20];
var x = new Array(); // 出力配列（解が入る配列）はメイン側で定義しておく
x = GaussJordan(a, b, n);
alert("nを与えたとき：x[0]=" + x[0] + ", x[1]=" + x[1] + ", x[2]=" + x[2]);
x = GaussJordan(a, b);
alert("nを省略したとき：x[0]=" + x[0] + ", x[1]=" + x[1] + ", x[2]=" + x[2]);
}
上の結果がalertされます。
データを入力してください↓
行・列のサイズn の配列 a[i][j] の逆行列を、ガウス=ジョルダンの掃出法で求めます。メイン側で定義した二次元配列 ｘ を与えれば、x に逆行列が入り、a は壊されません。ｘを省略すると a が逆行列になります。
n = 3; invMatrix(n, a, x);のとき invMatrix(n, a);のとき
a[0] = [ 1, 2. 1]; a[0} = [ 1, 2. 1] x[0] = [ -0.4, 0.6, 0.2] a[0] = [ -0.4, 0.6, 0.2]
a[1] = [ 2, 1. 0]; → a[1] = [ 2, 1. 0] x[1] = [ 0.8,-0.2,-0.4] a[1] = [ 0.8,-0.2,-0.4]
a[2] = [ 1, 1, 2]; a[2] = [ 1, 1, 2] x[2] = [ -0.2,-0.2, 0.6] a[2] = [ -0.2,-0.2, 0.6]
x = new Array();が必要か？ 必要 不要
トレース版を使うときは、
script： var ary = LPTrace("〇〇", m, n, a);
HTML： <div id="〇〇"></div>
とします。
シンプレックス法により線形計画法を解きます。それに加えてシャドウプライス（レジューズドコスト）や価格・定数項のレンジなどの感度分析もしています。参照：「シンプレックス法」、「シャドウ・プライスとレンジ」
変数
０ １ ２ ... ｎ
┌─┬─────────┐ ｍ，ｎ，ａは左図を参照してください。
０│ │ 目的関数 │
├─┼─────────┤ a[0][0] には "max"(最大化) "min"(最小化)が入ります。
１│ │ │
制 ２│定│ │ 定数項はすべて正数で、定数項≧～ の形式だけを扱います。
約 │数≧ │
式 ：│項│ │ 誤差の累積などは考慮していません。
ｍ│ │ │ 小規模な問題（ｍもｎも１０以下）を想定しています。
└─┴─────────┘
関数内部では、次のようなａを拡張した作業用配列を用いるので、ａの値は保持されます。
nx+1 n=nx+m
0 1 2 ↓ 3 4 5 ↓
+------+-------------+-------------------+-----
0 | | | | :←目的関数
+------+-------------+-------------------+ :
1 | | | | :
| 定 | 基幹部 | スラック部 | m+1
2 | 数 | | | :
| 項 | | | :
m = 3 | | | | :
+------+-------------+-------------------+-----
|----- nx ----|-------- m --------|
|-------------- nx + m + 1 -------------|
解は２次元配列になります。
メイン側で、 ary = LP(m, n, a); とすると、
ary の内容は次のようになります。
ary[3][1] = 72, ary[2][5] = 3
↓j k→ 0 1 2 3 4 5 0 変数名 目的関数 解の値 レジューズド レンジ レンジ 定数項 コスト 下限 上限 1 Ｘ1 3 16 0 2 8 2 Ｘ2 2 8 0 0.75 3 3 λ1 72 0 0.333 60 96 4 λ2 48 0 0.833 36 52.364 5 λ3 48 8 0 -0.5 0.455
変数名は、基幹変数（与えた変数）はＸj、制約式（スラック変数）はλiになります。
基底に入らないＸ、基底に入ったλのレンジは計算していません。
供給地ｉの供給量a[i][0] (i=1～m), 需要地ｊの需要量a[0][j] (j=1～n) と、ｉ→ｊの単位輸送費用a[i][j]を与えて、総輸送費用を最小化するように、ｉ→ｊの最適輸送量を計算します。
北西隅の方法により初期解を得て、Ｕ－Ｖ法により解の改善を行い、最適解を求めます。
すべての入力値は正の整数で、Σa[i][0] = Σa[0][j] でなければなりません。
入力データ
供給地数(m) = 3， 需要地数(n) = 4
a[i][j] ｊ→ ０ １ ２ ３ ４=ｎ
需要地1 需要地2 需要地3 需要地4
↓ｉ 12 20 16 18 ←需要量
１ 供給地1 18 7 3 2 10
２ 供給地2 22 9 3 6 8
３ 供給地3 26 8 7 8 6
↑ｍ ↑供給量
関数呼び出し
var ary = TLP(m, n, a);
トレース版を使うときは
script： var ary = TLPTrace("〇〇", m, n, a);
HTML： <div id="〇〇"></div>
戻り値（aryの内容)は、ary[i][j]の２次元配列になります。
i→ 0 1 2
↓ 供給地 需要地 輸送量
0 1 1 2
1 1 3 16
2 2 1 2
3 2 2 20
4 3 1 8
5 3 4 18
↑
ary[0].length = 6
総費用 = Σ(a[ary[0][j]][ary[1][j]]*ary[2][j]) = 296
基本的には、a, b, zは複素数である。a(1+2i) はa[0]=1; a[1]=2; b(3+4i) はb[0]=3; b[1]=4; のように、実数部を
[0]、虚数部を[1] とする配列として表現する。
単項演算：z = -a は z = complex("-", a); と表記する。zは(-1-2i) であるから、z[0]=-1, z[1]=-2 となる。
二項演算：z = a+b は、z = complex(a, "+", b); と表記する。z[0]=4, z[1]=6 となる。
この "-" や "+" などの演算子を op と表記する。
プログラムの記述
var a = [];
var z; // ★でｚの型が決まるので、z = []; とはしない。
a[0] = 1;
a[1] = 2;
z = complex("-", a); // ★
// z[0] = -1; z[1] = -2 になっている。
以下の一般形式での変数の型
a,b 入力 複素数あるいは実数スカラー
d 入力、出力 >複素数配列
z 出力 複素数
x 出力 実数スカラー
n 入力 正整数
基本統計量を得る。
リスト入力の場合 var ary = basicStat(40, 20, 50, 10, 30);
配列入力の場合 var x = [40, 20, 50, 10, 30];
var ary = basicStatArray(x); 配列ｘの全要素が対象
var ary = basicStatArray(x, 4); 先頭の４要素が対象（30 は対象外）
ary は統計量名称が要素とした連想配列となる。
ary["標本数"} = 5 n
ary["合計"} = 150 Σx
ary["平方和"} = 5500 Σx2
ary["平均"} = 30 μ = Σx／n
ary["偏差平方和"} = 1000 Σ(x-μ)2 = Σx2 - (Σx)2/n<
ary["不偏分散"} = 250 σ2 = 偏差平方和/(n-1)
ary["標準偏差"} = 15.81 σ = √不偏分散
ary["変動係数"} = 0.527 標準偏差/μn
ary["最小値"} = 10
ary["最大値"} = 50
ary["範囲"} = 40 最大値 - 最小値
ary["中央値"} = 30 メディアン
参照：「平均、分散、標準偏差」
基本統計量を得る。
リスト入力の場合 var ary = basicStat2(26,54, 25,62, 23,51, 27,58, 28,63,
25,65, 22,59, 26,63, 25,65, 23,60);
配列入力の場合 var x = [26, 25, 23, 27, 28, 25, 22, 26, 25, 23];
var y = [54, 62, 51, 58, 63, 65, 59, 63, 65, 60];
var ary = basicStatArray(x); 配列ｘの全要素（１０組）が対象
var ary = basicStatArray(x, 8); 先頭の８要素が対象（(25,65) (23,50) は対象外）
ary は統計量名称が要素とした連想配列となる。
ary["標本数"] = 10 n
ary["平均Ｘ"] = 25 μx = Σ(x)/n
ary["平均Ｙ"] = 60 μy = Σ(y)/n
ary["平方和Ｘ"] = 6282 Σ(x2)
ary["平方和Ｙ"] = 36194 Σ(y2)
ary["積和"] = 15023 Σ(x・y)
ary["偏差平方和Ｘ"] = 32 Sxx = Σ(x-μx)2
ary["偏差平方和Ｙ"] = 194 Syy = Σ(y-μy)2
ary["偏差積和"] = 23 Sxy = Σ(x-μx)(y-μy)
ary["共分散"] = 2.556 Sxy/(n-1)
ary["回帰係数Ｘ"] = 0.718 Bx = Sxy/Sxx y = a + bx のｂ
ary["回帰係数Ｙ"] = 0.119 By = Sxy/Syy x = a + by のｂ
ary["寄与率"] = 0.085 r2 = Bx・By
ary["相関係数"] = 0.291 r = √(r2)
Ｘ、Ｙは大文字全角
平均μ、標準偏差σの正規分布において、値ｘの確率密度関数の値ｐを戻します。
p = 1/√(２πσ2) * e(-(x-μ)2/2σ2)
μ、σ2を省略すると、μ＝０、σ2＝１の標準正規分布になります。
平均μ、分散σ2の正規分布において、－∞～ｘの確率密度関数の値ｓを戻します。μ、σ2を省略すると、μ＝０、σ2＝１の標準正規分布 sNormalDistStd() になります。
上の一般形の数値での意味：男性成人の平均身長は170cmで標準偏差が10cmであるとする。160cm以下の割合は？ 解：15.9%
sNormalDist()は、確率密度関数 pNormalDist() を数値積分 sFx() で計算してもよいのですが、直接計算する近似式が知られています。ここではHastingsの式を使いました。 （引用先：http://qiita.com/gigamori/items/e17e6f9faffb78822c56 by Yanai Takamichi）
正規分布（μ、σ）での xmin～xmax の累積確率 s を戻します。
例えば、「身長が平均（μ）170cm、標準偏差（σ）10cmの正規分布に従うとき、身長が160cm～180cmの割合は？ → 答 = 68.3%
正規分布の累積確率関数 sNormalDist() の逆関数です。累積確率がｐになる点ｘを戻します。
例えば、「身長が平均（μ）170cm、標準偏差（σ）10cmの正規分布に従うとき、確率（ｐ）90%の人はｘcm以下である」かを求めるとｘ=183cmだとなります。
標準形 invNormalDistStd() は、μ＝０、σ＝１の場合です。 sNormalDist() でμとσを省略した場合も同じです。
p=0.01 → x=-2.32, p=0.05 → x=-1.65, p=0.50 → x=0, p=0.55 → x=1.65, p=0.99 → x=2.32
となります（p<0.5 のときは x が負になります）。
invNormalDistStd() での逆関数の近似式は、「https://potatodigger.wordpress.com/2013/06/30/正規分布の確率密度関数pdf、累積分布関数cdf、逆累/ by potatodigger」を参考にさせていただきました。
| 関数一般形式 | 機能 | 計算式 | ソースリスト |
| p = pExpDist(t, λ) | t=t における確率密度 | p = λ*e(-λ*t) | |
| s = sExpDist(t, λ) | 区間 t=0～t の累積確率 | s = 1-e(-λ*t) | |
| s = sExpDistMinmax(tmin, tmax, λ) | 区間 t=tmin～tmax の累積確率 | s = e(-λ*tmin)-e(-λ*tmax) | |
| t = invExpDist(p, λ) | 累積確率が p となる t の値 | t = -(1/λ)log(1-p) |
指数分布と一様分布の中間的な分布。ｋ＝１なら指数分布、ｋ＝∞なら一様分布
| 関数一般形式 | 機能 | 計算式 | ソースリスト |
| p = pErlangDist(t, λ, k) | t=t における確率密度 | p = (λ*t)k-1/(k-1)! * λe(-λ*t) | |
| s = sErlangDist(t, λ, k) | 区間 t=0～t の累積確率 | pErlangDistの単純累計 |
λは正の実数、ｋは自然数、ｔは本来正の実数だが、ここでは自然数に限定する
密度確率（ｔ＝ｔでの確率）
累積確率（０～ｔの累積）
λとｋを与えて、ｔ＝０～tmaxでの確率密度と累積密度の表を作成
一定（微小）時間内にある事象が発生する回数が平均λであるとき、発生回数がｘとなる確率ｐは、
p(x) = (λx/x!)*e-λ ← pPoissonDist(x,λ)
のポアソン分布に従う（参照：「ポアソン分布と指数分布」
累積確率は、単純に s(x) = s(x-1) + p(x) で求めた。
逆関数は、x=0, 1, 2, ... として累積確率 s(x) を計算し、s(x) ≦ p である最大のｘを戻している。すなわち「確率ｐでｘ回以下だといえる」こととした。
λは正の実数、pは０≦ｐ≦１、ｘは０以上の整数でなければならない。
発生確率がｐの事象が、ｎ回の試行でｘ回発生する確率。
確率密度： ｘ回発生する確率 Ｐ(x) = nＣx・px(1-p)n-x = pBinDist(x, n, p)
累積密度： ０～ｘ回発生する確率 Ｓ(≦x) = ΣＰ(x) = sBinDist(x, n, p)
nＣx = n!/x!(n-x)! = nCx(n, x) 二項係数（ｎ個からｘ個を取り出す場合の数（場合の数））
参照：「順列・組合せ」、
「二項分布」、
「２項分布の計算プログラム」
ここの例では、ｘを与えず、ｘ＝０～ｎの表を掲げます。
二項分布の計算は計算量が大きく、大きい数と小さい数の計算が必要になるので誤差が生じます。計算を容易にするために、
・np > 5 かつ n(1-p) > 5 ならば、平均μ=np、標準偏差σ=√np(1-p) の正規分布で近似できる。
確率密度：区間 x-0.5～x+0.5 の正規分布累積確率 sNormalDistMinmax(x-0.5, x+0.5, μ, σ)
累積密度：≦x+0.5 の累積確率 sNormalDist(x+0.5, μ, σ)
・n>50, np≦5 ならば、平均λ=np のポアソン分布で近似できる。
などの方法がありますが、ここでは用いていません。
負の二項分布とは、ある事象が起こる確率がｐのとき、ｎ回試行したときに、その事象がはじめてｘ回になる確率Ｐ(n)の確率分布です。
Ｐ(n) = n-1Ｃx-1・ｐx(1-p)n-x
で計算できます（参照：「二項分布」）。
x=1 のときは、幾何分布 pGeoDist と一致します
発生確率ｐの事象がはじめて発生したときの試行回数ｎの確率分布（負の二項分布で、x＝１とおくと、幾何分布に一致）
確率密度： p = pGeoDist(n, p) = (1-p)n-1*ｐ
累積確率： s = sGeoDist(n, p) = 1-(1-p)n
参照：「二項分布」
確率密度：袋の中にＡがａ個、Ｂがｂ個入っている。α＋β個を取り出したとき、Ａがα個、Ｂがβ個である確率
p = aＣα・bＣβ／a+bＣα+β
で計算できる（参照：「二項分布」）。
確率密度：αが０～α個である確率。「α以上である確率」は sHyperGeoDist(a, b, α-1, β+1) になる。
３点 (xmin,0), (xmod, h), (xmax,0) を結ぶ三角形となる分布です。
ここでは、いろいろな点の値を計算して表にしています。
ｘを整数値にしたいときがあります。連続量の場合は、xmin および xmaxでは発生確率が０になりますが、整数の場合はその点でも０ではない値であるのが通常です。そのときは、やや乱暴ですが、xmin-0.5～xmax+0.5とします。そして累積確率では、x+0.5 の点を累積確率とし、逆関数のときは得られた値を四捨五入ば、それらしい値が得られます。
下の表では、連続のときは、得たｘを四捨五入して表示しているので、最小値と最大値の割合は他の半分になっています。
整数のときは、区間を xmin～xmin+1として、乱数を発生させ、得たｘを切り捨てているので、同じ値になっています。
0～1の値を持つ乱数から、正規分布の累積確率の逆関数を用いて、値ｘを得ます。
整数値のときは、単に四捨五入しているだけです。区間指定ありの場合は、区間外の値になったときは。区間内になるまで繰返し乱数を発生させます。
xRandomNormal12はxRandomNormalと全く同じ機能ですが、考え方が異なります。
中心極限定理により、多数の一様分布に従う標本の和は正規分布になります。Math.random()は０～１の一様分布です。それを１２個を合計した値は、μ＝６、σ＝１の正規分布に近づきます。xRandomNormal(μ, σ)が内部でinvNormalDistを用いているのに対して、xRandomNormal12は非常に簡単なアルゴリズムで実現しています。精度は、xRandomNormalのほうが良いようです。
単位時間中にある事象が発生する平均回数をλとするとき、その事象の発生間隔がx単位時間である確率は指数分布に従う。
指数分布累積確率の逆関数は、x = -1/λ*log(1-p) であるから、指数分布乱数は、x = -1/λ*log(1-Math.random()) で得られる。ｘの平均値は 1/λ になる。
なお、事象発生は稀であるとの前提から、λの値が小さくなるよう単位時間を設定することが望ましい。
参照：「ポアソン分布と指数分布」
ここでの出力について
xmodは最頻値、平均は (xmin+xmax+xmod)/3 になる。
三角分布の累積確率の逆関数を利用してもよいが、二つの一様分布乱数、u, v (u>vとする）として、
(1-c)*v + c*u c = (xmod-xmin)/(xmax-xmin)
とすることにより、xmin=0, xmax=1, xmod=c の三角分布乱数が得られる。このほうが簡単なので、ここではこれを採用した。