Module: Numo::GSL::Multifit

Defined in:
ext/numo/gsl/multifit/gsl_multifit.c

Defined Under Namespace

Classes: LinearWorkspace

Class Method Summary collapse

Class Method Details

.linear(X, y) ⇒ GSL::Multifit::LinearResult

This function computes the best-fit parameters c of the model y = X c for the observations y and the matrix of predictor variables X, using the preallocated workspace provided in work. The p-by-p variance-covariance matrix of the model parameters cov is set to \sigma^2 (X^T X)^-1, where \sigma is the standard deviation of the fit residuals. The sum of squares of the residuals from the best-fit, \chi^2, is returned in chisq. If the coefficient of determination is desired, it can be computed from the expression R^2 = 1 - \chi^2 / TSS, where the total sum of squares (TSS) of the observations y may be computed from gsl_stats_tss.

The best-fit is found by singular value decomposition of the matrix X using the modified Golub-Reinsch SVD algorithm, with column scaling to improve the accuracy of the singular values. Any components which have zero singular value (to machine precision) are discarded from the fit.

Parameters:

  • X (DFloat) —

    (input matrix) predictor variables

  • y (DFloat) —

    (input vector) observation

Returns:

  • (GSL::Multifit::LinearResult) —

    result Struct with members: c, cov, chisq.



129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
 
# File 'ext/numo/gsl/multifit/gsl_multifit.c', line 129

static VALUE
multifit_linear_workspace_s_linear(VALUE mod, VALUE v1, VALUE v2)
{
    size_t c_shape[1], cov_shape[2];
    ndfunc_arg_in_t ain[2] = {{cDF,2},{cDF,1}};
    ndfunc_arg_out_t aout[3] = {{cDF,1,c_shape},{cDF,2,cov_shape},{cDF,0}};
    ndfunc_t ndf = { iter_multifit_linear_workspace_s_linear, NO_LOOP|NDF_EXTRACT,
                     2, 3, ain, aout };
    VALUE vws, r, result;
    narray_t *X, *y;
    size_t n, p;
    gsl_multifit_linear_workspace *ws;

    GetNArray(v1,X);
    GetNArray(v2,y);
    CHECK_GE_2D(X);
    CHECK_GE_1D(y);
    n = MAT_SIZE1(X);
    p = MAT_SIZE2(X);
    CHECK_SIZE_EQ(n,VEC_SIZE(y),"y size does not match X column size");

    c_shape[0] = cov_shape[0] = cov_shape[1] = p;

    ws = gsl_multifit_linear_alloc(n,p);
    if (!ws)
        rb_raise(rb_eNoMemError,"fail to allocate struct");
    vws = TypedData_Wrap_Struct(cLinearWorkspace, &multifit_linear_workspace_data_type, (void*)ws);

    r = na_ndloop3(&ndf, ws, 2, v1, v2);

    result = rb_class_new_instance(3, RARRAY_PTR(r), cLinearResult);

    RB_GC_GUARD(vws);
    RB_GC_GUARD(r);
    return result;
}

.linear_est(x, c, cov) ⇒ [DFloat,DFloat]

This function uses the best-fit multilinear regression coefficients c and their covariance matrix cov to compute the fitted function value y and its standard deviation y_err for the model y = x.c at the point x.

Parameters:

  • x (DFloat) —

    (input vector) observations

  • c (DFloat) —

    (input vector) solusions

  • cov (DFloat) —

    (input matrix) covariance

Returns:

  • ([DFloat,DFloat]) —

    array of (y, y_err).



282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
 
# File 'ext/numo/gsl/multifit/gsl_multifit.c', line 282

static VALUE
multifit_s_linear_est(VALUE mod, VALUE v1, VALUE v2, VALUE v3)
{
    ndfunc_arg_in_t ain[3] = {{cDF,1},{cDF,1},{cDF,2}};
    ndfunc_arg_out_t aout[2] = {{cDF,0},{cDF,0}};
    ndfunc_t ndf = { iter_multifit_s_linear_est, NO_LOOP|NDF_EXTRACT,
                     3, 2, ain, aout };
    narray_t *x, *c, *cov;

    GetNArray(v1,x);
    GetNArray(v2,c);
    GetNArray(v3,cov);
    CHECK_GE_1D(x);
    CHECK_GE_1D(c);
    CHECK_GE_2D(cov);
    CHECK_SIZE_EQ(VEC_SIZE(x),VEC_SIZE(c),"x size does not match c size");
    CHECK_SIZE_EQ(MAT_SIZE1(cov),MAT_SIZE2(cov),"cov is not square matrix");
    CHECK_SIZE_EQ(VEC_SIZE(c),MAT_SIZE1(cov),"c size does not match cov row size");
    return na_ndloop(&ndf, 3, v1, v2, v3);
}

.linear_residuals(X, y, c) ⇒ DFloat

This function computes the vector of residuals r = y - X c for the observations y, coefficients c and matrix of predictor variables X.

Parameters:

  • X (DFloat) —

    (input matrix) design matrix

  • y (DFloat) —

    (input vector) rhs vector

  • c (DFloat) —

    (input vector) fit coefficients

Returns:

  • (DFloat) —

    returns r



333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
 
# File 'ext/numo/gsl/multifit/gsl_multifit.c', line 333

static VALUE
multifit_s_linear_residuals(VALUE mod, VALUE v1, VALUE v2, VALUE v3)
{
    size_t shape[1];
    ndfunc_arg_in_t ain[3] = {{cDF,2},{cDF,1},{cDF,1}};
    ndfunc_arg_out_t aout[1] = {{cDF,1,shape}};
    ndfunc_t ndf = { iter_multifit_s_linear_residuals, NO_LOOP|NDF_EXTRACT,
                     3, 1, ain, aout };
    narray_t *x, *y, *c;

    GetNArray(v1,x);
    GetNArray(v2,y);
    GetNArray(v3,c);
    CHECK_GE_2D(x);
    CHECK_GE_1D(y);
    CHECK_GE_1D(c);
    CHECK_SIZE_EQ(MAT_SIZE1(x),VEC_SIZE(y),"y size does not match x row size");
    CHECK_SIZE_EQ(MAT_SIZE2(x),VEC_SIZE(c),"c size does not match x column size");
    shape[0] = VEC_SIZE(y);
    return na_ndloop(&ndf, 3, v1, v2, v3);
}

.wlinear(X, w, y) ⇒ GSL::Multifit::LinearResult

This function computes the best-fit parameters c of the weighted model y = X c for the observations y with weights w and the matrix of predictor variables X, using the preallocated workspace provided in work. The p-by-p covariance matrix of the model parameters cov is computed as (X^T W X)^-1. The weighted sum of squares of the residuals from the best-fit, \chi^2, is returned in chisq. If the coefficient of determination is desired, it can be computed from the expression R^2 = 1 - \chi^2 / WTSS, where the weighted total sum of squares (WTSS) of the observations y may be computed from gsl_stats_wtss.

Parameters:

  • X (DFloat) —

    (input matrix) predictor variables

  • w (DFloat) —

    (input vector) observations

  • y (DFloat) —

    (input vector) weights

Returns:

  • (GSL::Multifit::LinearResult) —

    result Struct with members: c, cov, chisq.



208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
 
# File 'ext/numo/gsl/multifit/gsl_multifit.c', line 208

static VALUE
multifit_linear_workspace_s_wlinear(VALUE mod, VALUE v1, VALUE v2, VALUE v3)
{
    size_t c_shape[1], cov_shape[2];
    ndfunc_arg_in_t ain[3] = {{cDF,2},{cDF,1},{cDF,1}};
    ndfunc_arg_out_t aout[3] = {{cDF,1,c_shape},{cDF,2,cov_shape},{cDF,0}};
    ndfunc_t ndf = { iter_multifit_linear_workspace_s_wlinear, NO_LOOP|NDF_EXTRACT,
                     3, 3, ain, aout };
    VALUE vws, r, result;
    narray_t *X, *y, *w;
    size_t n, p;
    gsl_multifit_linear_workspace *ws;

    GetNArray(v1,X);
    GetNArray(v2,y);
    GetNArray(v3,w);
    CHECK_GE_2D(X);
    CHECK_GE_1D(y);
    CHECK_GE_1D(w);
    n = MAT_SIZE1(X);
    p = MAT_SIZE2(X);
    CHECK_SIZE_EQ(n,VEC_SIZE(y),"y size does not match X column size");
    CHECK_SIZE_EQ(n,VEC_SIZE(w),"w size does not match X column size");

    c_shape[0] = cov_shape[0] = cov_shape[1] = p;

    ws = gsl_multifit_linear_alloc(n,p);
    if (!ws)
        rb_raise(rb_eNoMemError,"fail to allocate struct");
    vws = TypedData_Wrap_Struct(cLinearWorkspace, &multifit_linear_workspace_data_type, (void*)ws);

    r = na_ndloop3(&ndf, ws, 3, v1, v2, v3);

    result = rb_class_new_instance(3, RARRAY_PTR(r), cLinearResult);

    RB_GC_GUARD(vws);
    RB_GC_GUARD(r);
    return result;
}