// Extend the Array class
Array.prototype.max = function () {
return Math.max.apply(null, this);
};
Array.prototype.min = function () {
return Math.min.apply(null, this);
};
Array.prototype.mean = function () {
var i, sum;
for (i = 0, sum = 0; i < this.length; i++)
{sum += this[i];}
return sum / this.length;
};
Array.prototype.pip = function (x, y) {
var i, j, c = false;
for (i = 0, j = this.length - 1; i < this.length; j = i++) {
if (((this[i][1] > y) != (this[j][1] > y)) &&
(x < (this[j][0] - this[i][0]) * (y - this[i][1]) / (this[j][1] - this[i][1]) + this[i][0])) {
c = !c;
}
}
return c;
}
var kriging = function () {
var kriging = {};
var createArrayWithValues = function (value, n) {
var array = [];
for (var i = 0; i < n; i++) {
array.push(value);
}
return array;
};
// Matrix algebra矩阵代数
kriging_matrix_diag = function (c, n) {
var Z = createArrayWithValues(0, n * n);
for (i = 0; i < n; i++) {Z[i * n + i] = c;}
return Z;
};
//将这个矩阵变为转置阵,也就是将元素颠倒顺序
kriging_matrix_transpose = function (X, n, m) {
var i, j, Z = Array(m * n);
for (i = 0; i < n; i++)
{
for (j = 0; j < m; j++)
{Z[j * n + i] = X[i * m + j];}
}
return Z;
};
//再次改变数值,把c给每一个二维元素赋值
kriging_matrix_scale = function (X, c, n, m) {
var i, j;
for (i = 0; i < n; i++)
{
for (j = 0; j < m; j++)
{X[i * m + j] *= c;}
}
};
kriging_matrix_add = function (X, Y, n, m) {
var i, j, Z = Array(n * m);
for (i = 0; i < n; i++)
{
for (j = 0; j < m; j++)
{Z[i * m + j] = X[i * m + j] + Y[i * m + j];}
}
return Z;
};
// Naive matrix multiplication
//简单的矩阵乘法,矩阵和矩阵的乘法
//也就是前一个矩阵中的行乘以后一个矩阵中的列
kriging_matrix_multiply = function (X, Y, n, m, p) {
var i, j, k, Z = Array(n * p);
for (i = 0; i < n; i++) {
for (j = 0; j < p; j++) {
Z[i * p + j] = 0;
for (k = 0; k < m; k++)
{Z[i * p + j] += X[i * m + k] * Y[k * p + j];}
}
}
return Z;
};
// Cholesky decomposition
//柯列斯基分解,这是一种将正定矩阵分解为上三角矩阵和下三角矩阵的方法,
//在优化矩阵计算的时候会用到的一种技巧
//也就是,下面左边为下三角,右边为上三角
kriging_matrix_chol = function (X, n) {
var i, j, k, sum, p = Array(n);
for (i = 0; i < n; i++) {p[i] = X[i * n + i];}
for (i = 0; i < n; i++) {
for (j = 0; j < i; j++){
p[i] -= X[i * n + j] * X[i * n + j];
}
if (p[i] <= 0) return false;
p[i] = Math.sqrt(p[i]);
for (j = i + 1; j < n; j++) {
for (k = 0; k < i; k++){
X[j * n + i] -= X[j * n + k] * X[i * n + k];
}
X[j * n + i] /= p[i];
}
}
for (i = 0; i < n; i++) X[i * n + i] = p[i];
return true;
};
// Inversion of cholesky decomposition
// 用斯基分解求矩阵的逆
kriging_matrix_chol2inv = function (X, n) {
var i, j, k, sum;
for (i = 0; i < n; i++) {
X[i * n + i] = 1 / X[i * n + i];
for (j = i + 1; j < n; j++) {
sum = 0;
for (k = i; k < j; k++){
sum -= X[j * n + k] * X[k * n + i];
}
X[j * n + i] = sum / X[j * n + j];
}
}
for (i = 0; i < n; i++){
for (j = i + 1; j < n; j++){
X[i * n + j] = 0;
}
}
for (i = 0; i < n; i++) {
X[i * n + i] *= X[i * n + i];
for (k = i + 1; k < n; k++){
X[i * n + i] += X[k * n + i] * X[k * n + i];
}
for (j = i + 1; j < n; j++){
for (k = j; k < n; k++){
X[i * n + j] += X[k * n + i] * X[k * n + j];
}
}
}
for (i = 0; i < n; i++){
for (j = 0; j < i; j++){
X[i * n + j] = X[j * n + i];
}
}
};
// Inversion via gauss-jordan elimination
//用高斯-约当消去法求逆,它的速度不是最快的,但是它非常稳定
//如果A是求解矩阵,那么求A的逆矩阵则为
//用A矩阵右边乘以单位矩阵I(与A同行同列值为1的单位矩阵)
//公式为A*I=I*B,(等号右边要同时变化),也就是一个矩阵右边乘以单位矩阵化为,
//左边单位矩阵乘以B,则B就是A矩阵的逆
kriging_matrix_solve = function (X, n) {
var m = n;
var b = Array(n * n);
var indxc = Array(n);
var indxr = Array(n);
var ipiv = Array(n);
var i, icol, irow, j, k, l, ll;
var big, dum, pivinv, temp;
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++) {
if (i == j) b[i * n + j] = 1;
else b[i * n + j] = 0;
}
}
for (j = 0; j < n; j++){ ipiv[j] = 0;}
for (i = 0; i < n; i++) {
big = 0;
for (j = 0; j < n; j++) {
if (ipiv[j] != 1) {
for (k = 0; k < n; k++) {
if (ipiv[k] == 0) {
if (Math.abs(X[j * n + k]) >= big) {
big = Math.abs(X[j * n + k]);
irow = j;
icol = k;
}
}
}
}
}
++(ipiv[icol]);
if (irow != icol) {
for (l = 0; l < n; l++) {
temp = X[irow * n + l];
X[irow * n + l] = X[icol * n + l];
X[icol * n + l] = temp;
}
for (l = 0; l < m; l++) {
temp = b[irow * n + l];
b[irow * n + l] = b[icol * n + l];
b[icol * n + l] = temp;
}
}
indxr[i] = irow;
indxc[i] = icol;
if (X[icol * n + icol] == 0) return false; // Singular
pivinv = 1 / X[icol * n + icol];
X[icol * n + icol] = 1;
for (l = 0; l < n; l++){ X[icol * n + l] *= pivinv;}
for (l = 0; l < m; l++) {b[icol * n + l] *= pivinv;}
for (ll = 0; ll < n; ll++) {
if (ll != icol) {
dum = X[ll * n + icol];
X[ll * n + icol] = 0;
for (l = 0; l < n; l++) X[ll * n + l] -= X[icol * n + l] * dum;
for (l = 0; l < m; l++) b[ll * n + l] -= b[icol * n + l] * dum;
}
}
}
for (l = (n - 1); l >= 0; l--)
{
if (indxr[l] != indxc[l]) {
for (k = 0; k < n; k++) {
temp = X[k * n + indxr[l]];
X[k * n + indxr[l]] = X[k * n + indxc[l]];
X[k * n + indxc[l]] = temp;
}
}
}
return true;
}
// Variogram models
//变差函数模型
//变差函数高斯
kriging_variogram_gaussian = function (h, nugget, range, sill, A) {
return nugget + ((sill - nugget) / range) *
(1.0 - Math.exp(-(1.0 / A) * Math.pow(h / range, 2)));
};
// 变差函数指数
kriging_variogram_exponential = function (h, nugget, range, sill, A) {
return nugget + ((sill - nugget) / range) *
(1.0 - Math.exp(-(1.0 / A) * (h / range)));
};
// 变差函数的球形
kriging_variogram_spherical = function (h, nugget, range, sill, A) {
if (h > range) return nugget + (sill - nugget) / range;
return nugget + ((sill - nugget) / range) *
(1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
};
// Train using gaussian processes with bayesian priors
//训练使用高斯过程与贝叶斯先验
//kriging.train(t, x, y, model, sigma2, alpha):
//使用gaussian、exponential或spherical模型对数据集进行训练,返回的是一个variogram对象;
kriging.train = function (t, x, y, model, sigma2, alpha) {
var variogram = {
t: t,
x: x,
y: y,
nugget: 0.0,
range: 0.0,
sill: 0.0,
A: 1 / 3,
n: 0
};
switch (model) {
case "gaussian":
variogram.model = kriging_variogram_gaussian;
break;
case "exponential":
variogram.model = kriging_variogram_exponential;
break;
case "spherical":
variogram.model = kriging_variogram_spherical;
break;
};
// Lag distance/semivariance
// 滞后距离/半方差
var i, j, k, l, n = t.length;
var distance = Array((n * n - n) / 2);
for (i = 0, k = 0; i < n; i++)
{
for (j = 0; j < i; j++, k++) {
distance[k] = Array(2);
distance[k][0] = Math.pow(
Math.pow(x[i] - x[j], 2) +
Math.pow(y[i] - y[j], 2), 0.5);
distance[k][1] = Math.abs(t[i] - t[j]);
}
}
distance.sort(function (a, b) { return a[0] - b[0]; });
variogram.range = distance[(n * n - n) / 2 - 1][0];
// Bin lag distance 本滞后距离
var lags = ((n * n - n) / 2) > 30 ? 30 : (n * n - n) / 2;
var tolerance = variogram.range / lags;
var lag = createArrayWithValue
- 1
- 2
- 3
- 4
- 5
- 6
前往页