Skip to content

[Backport 6.1] Calculate confidence band for -Er and -Eo #3655

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jul 13, 2020
Merged

[Backport 6.1] Calculate confidence band for -Er and -Eo #3655

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 44 additions & 10 deletions src/gmtregress.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,8 @@ enum GMT_enum_regress {
GMTREGRESS_YMEAN = 7,
GMTREGRESS_R = 8,
GMTREGRESS_CORR = 9,
GMTREGRESS_NPAR = 10,
GMTREGRESS_MISFTY = 10,
GMTREGRESS_NPAR = 11,
GMTREGRESS_NPAR_MAIN = 4,
GMTREGRESS_OUTPUT_GOOD = 1,
GMTREGRESS_OUTPUT_BAD = 2};
Expand Down Expand Up @@ -683,7 +684,7 @@ GMT_LOCAL double gmtregress_LSy_regress1D (struct GMT_CTRL *GMT, double *x, doub
par[GMTREGRESS_SIGIC] = sqrt (S_xx * D);
for (k = 0; k < n; k++) /* Here we recycle Q to hold y-residual e */
Q[k] = y[k] - gmtregress_model (x[k], par);
par[GMTREGRESS_MISFT] = gmtregress_L2_misfit (GMT, Q, W, n, GMTREGRESS_Y, 0.0);
par[GMTREGRESS_MISFT] = par[GMTREGRESS_MISFTY] = gmtregress_L2_misfit (GMT, Q, W, n, GMTREGRESS_Y, 0.0);
par[GMTREGRESS_ANGLE] = atand (par[GMTREGRESS_SLOPE]);
scale = gmtregress_L2_scale (GMT, NULL, W, n, par);
gmt_M_free (GMT, Q);
Expand All @@ -693,8 +694,32 @@ GMT_LOCAL double gmtregress_LSy_regress1D (struct GMT_CTRL *GMT, double *x, doub
return (scale);
}

GMT_LOCAL void gmtregress_yorkRMA_error (struct GMT_CTRL *GMT, double *U, double *V, uint64_t n, double sx, double sy, double *par) {
uint64_t k;
double *u = gmt_M_memory (GMT, NULL, n, double), *x = gmt_M_memory (GMT, NULL, n, double);
double v, sum_u2, mean_x;
gmt_M_unused (GMT);

/* From York et al [2004] for RMA case when w(Xi) = 1/sx^2 and w(Yi) = 1/sy2 */
sx *= sx; sy *= sy; /* Get variances */
v = sy + par[GMTREGRESS_SLOPE] * par[GMTREGRESS_SLOPE] * sx;
for (k = 0; k < n; k++)
x[k] = par[GMTREGRESS_XMEAN] + (U[k] * sy + par[GMTREGRESS_SLOPE] * V[k] * sx) / v;
mean_x = gmtregress_gmt_sum (x, n) / n; /* Get sum of x divided by n */
for (k = 0; k < n; k++) /* compute u */
u[k] = x[k] - mean_x;
gmtregress_eval_product (u, u, x, n); /* Compute x[i] = u[i] * u[i] */
sum_u2 = gmtregress_gmt_sum (x, n); /* Get sum of u^2 */
par[GMTREGRESS_SIGSL] = 1 / sum_u2;
par[GMTREGRESS_SIGIC] = (1 / n + mean_x * mean_x / sum_u2);
par[GMTREGRESS_SIGSL] = sqrt (par[GMTREGRESS_SIGSL]);
par[GMTREGRESS_SIGIC] = sqrt (par[GMTREGRESS_SIGIC]);
gmt_M_free (GMT, x);
gmt_M_free (GMT, u);
}

GMT_LOCAL double gmtregress_LSxy_regress1D_basic (struct GMT_CTRL *GMT, double *x, double *y, uint64_t n, double *par) {
/* Basic LS xy orthogonal regression, with no data errors. See York [1966] */
/* Basic LS xy orthogonal regression, with no data errors. See York [1966] and York et al [2004] */
uint64_t k;
unsigned int p;
double *u = gmt_M_memory (GMT, NULL, n, double), *v = gmt_M_memory (GMT, NULL, n, double);
Expand All @@ -712,7 +737,6 @@ GMT_LOCAL double gmtregress_LSxy_regress1D_basic (struct GMT_CTRL *GMT, double *
sum_v2 = gmtregress_gmt_sum (Q, n); /* Get sum of v*v */
gmtregress_eval_product (u, v, Q, n); /* Compute Q[i] = u[i] * v[i] */
sum_uv = gmtregress_gmt_sum (Q, n); /* Get sum of u*v */
gmt_M_free (GMT, u); gmt_M_free (GMT, v); /* Done with these arrays */
part1 = sum_v2 - sum_u2;
part2 = sqrt (pow (sum_u2 - sum_v2, 2.0) + 4.0 * sum_uv * sum_uv);
b[0] = (part1 + part2) / (2.0 * sum_uv);
Expand All @@ -727,39 +751,49 @@ GMT_LOCAL double gmtregress_LSxy_regress1D_basic (struct GMT_CTRL *GMT, double *
p = (E[0] < E[1]) ? 0 : 1; /* Determine the solution with the smallest misfit and copy to par array: */
par[GMTREGRESS_SLOPE] = b[p];
par[GMTREGRESS_ICEPT] = a[p];
par[GMTREGRESS_SIGSL] = par[GMTREGRESS_SLOPE] * sqrt ((1.0 - r * r) / n) / r;
par[GMTREGRESS_SIGIC] = sqrt (pow (sig_y - sig_x * par[GMTREGRESS_SLOPE], 2.0) / n + (1.0 - r) * par[GMTREGRESS_SLOPE] * (2.0 * sig_x * sig_y + (mean_x * par[GMTREGRESS_SLOPE] * (1.0 + r) / (r * r))));
gmtregress_yorkRMA_error (GMT, u, v, n, sig_x, sig_y, par);
par[GMTREGRESS_MISFT] = E[p];
/* Compute regular y LS misfit to use with confidence band */
for (k = 0; k < n; k++) Q[k] = y[k] - b[p] * x[k] - a[p];
par[GMTREGRESS_MISFTY] = gmtregress_L2_misfit (GMT, Q, W, n, GMTREGRESS_Y, 0.0);
par[GMTREGRESS_ANGLE] = atand (par[GMTREGRESS_SLOPE]);
par[GMTREGRESS_XMEAN] = mean_x;
par[GMTREGRESS_YMEAN] = mean_y;
scale = gmtregress_L2_scale (GMT, NULL, W, n, par);
gmt_M_free (GMT, Q);
gmt_M_free (GMT, W);
gmt_M_free (GMT, u);
gmt_M_free (GMT, v);

return (scale);
}

GMT_LOCAL double gmtregress_LSRMA_regress1D (struct GMT_CTRL *GMT, double *x, double *y, double *w[], uint64_t n, double *par) {
/* Basic LS RMA orthogonal regression with no weights [Reference?] */
/* Basic LS RMA orthogonal regression with no weights [York et al, 2004] */
uint64_t k;
double sx, sy, scale, r;
double *U = gmt_M_memory (GMT, NULL, n, double), *V = gmt_M_memory (GMT, NULL, n, double), *W = gmt_M_memory (GMT, NULL, n, double);
double sx, sy, scale, r, sum_u2;
double *U = gmt_M_memory (GMT, NULL, n, double), *V = gmt_M_memory (GMT, NULL, n, double), *Q = gmt_M_memory (GMT, NULL, n, double), *W = gmt_M_memory (GMT, NULL, n, double);
gmt_M_memset (par, GMTREGRESS_NPAR, double);
(void)gmtregress_demeaning (GMT, x, y, w, n, par, U, V, W, NULL, NULL);
r = gmt_corrcoeff (GMT, U, V, n, 1);
sx = gmt_std_weighted (GMT, U, w[GMT_X], 0.0, n);
sy = gmt_std_weighted (GMT, V, w[GMT_Y], 0.0, n);
gmtregress_eval_product (U, U, Q, n); /* Compute Q[i] = u[i] * u[i] */
sum_u2 = gmtregress_gmt_sum (Q, n); /* Get sum of u*u */
par[GMTREGRESS_SLOPE] = sy / sx;
if (r < 0.0) par[GMTREGRESS_SLOPE] = -par[GMTREGRESS_SLOPE]; /* Negative correlation means negative slope */
par[GMTREGRESS_ICEPT] = par[GMTREGRESS_YMEAN] - par[GMTREGRESS_SLOPE] * par[GMTREGRESS_XMEAN];
par[GMTREGRESS_ANGLE] = atand (par[GMTREGRESS_SLOPE]);
gmtregress_yorkRMA_error (GMT, U, V, n, sx, sy, par);
for (k = 0; k < n; k++) /* Here we recycle U as y-residual e */
U[k] = y[k] - gmtregress_model (x[k], par);
/* Report RMA misfit but use L2 y-misfit in calculations for confidence band */
par[GMTREGRESS_MISFT] = gmtregress_L2_misfit (GMT, U, W, n, GMTREGRESS_RMA, par[GMTREGRESS_SLOPE]);
par[GMTREGRESS_MISFTY] = gmtregress_L2_misfit (GMT, U, W, n, GMTREGRESS_Y, par[GMTREGRESS_SLOPE]);
scale = gmtregress_L2_scale (GMT, NULL, W, n, par);
gmt_M_free (GMT, U);
gmt_M_free (GMT, V);
gmt_M_free (GMT, Q);
gmt_M_free (GMT, W);
return (scale);
}
Expand Down Expand Up @@ -1336,7 +1370,7 @@ EXTERN_MSC int GMT_gmtregress (void *V_API, int mode, void *args) {
out[col] = S->data[GMT_Y][row] - gmtregress_model (x[row], par);
break;
case 'c': /* Model confidence limit (add slope and intercept uncertainties in quadrature since uncorrelated) */
out[col] = t_scale * sqrt (par[GMTREGRESS_MISFT]) * hypot (par[GMTREGRESS_SIGIC], par[GMTREGRESS_SIGSL] * fabs (x[row] - par[GMTREGRESS_XMEAN]));
out[col] = t_scale * sqrt (par[GMTREGRESS_MISFTY]) * hypot (par[GMTREGRESS_SIGIC], par[GMTREGRESS_SIGSL] * fabs (x[row] - par[GMTREGRESS_XMEAN]));
break;
case 'z': /* Standardized residuals (z-scores) */
out[col] = z_score[row];
Expand Down