Robust regression, compute regression with user-supplied functions and weights

はじめる前に
新機能一覧
Maple ワークシートの作成
Mapleワークシートを共有
Maple ウィンドウのカスタマイズ
Maple システムのカスタマイズ
コネクティビティ
Mathematics
数学
物理
プログラミング
グラフィックス
Student パッケージ
Science and Engineering
科学・エンジニアリング
アプリケーションと例題
リファレンス
システム
MapleSim
MapleSim ツールボックス
MapleSim ヘルプ
Tasks
Toolboxes
- BlockImporter
- Global Optimization
- Grid Computing
- Maple-NAG Connector
  - General information
  - Computational Routines
    - a02 (Complex Arithmetic)
    - c02 (Zeros of Polynomials)
    - c05 (Roots of One or More Transcendental Equations)
    - c06 (Fourier Transforms)
    - d01 (Quadrature)
    - d02 (Ordinary Differential Equations)
    - d03 (Partial Differential Equations)
    - d06 (Mesh Generation)
    - e01 (Interpolation)
    - e02 (Curve and Surface Fitting)
    - e04 (Minimizing or Maximizing a Function)
    - f01 (Matrix Factorizations)
    - f02 (Eigenvalues and Eigenvectors)
    - f03 (Determinants)
    - f04 (Simultaneous Linear Equations)
    - f06 (Linear Algebra Support Functions)
    - f07 (Linear Equations (LAPACK))
    - f08 (Least-squares and Eigenvalue Problems (LAPACK))
    - f11 (Large Scale Linear Systems)
    - f16 (NAG Interface to BLAS)
    - g01 (Simple Calculations on Statistical Data)
    - g02 (Correlation and Regression Analysis)
      - g02 Introduction
      - g02brc (nag_ken_spe_corr_coeff)
      - g02btc (nag_sum_sqs_update)
      - g02buc (nag_sum_sqs)
      - g02bwc (nag_cov_to_corr)
      - g02bxc (nag_corr_cov)
      - g02byc (nag_partial_corr)
      - g02cac (nag_simple_linear_regression)
      - g02cbc (nag_regress_confid_interval)
      - g02dac (nag_regsn_mult_linear)
      - g02dcc (nag_regsn_mult_linear_addrem_obs)
      - g02ddc (nag_regsn_mult_linear_upd_model)
      - g02dec (nag_regsn_mult_linear_add_var)
      - g02dfc (nag_regsn_mult_linear_delete_var)
      - g02dgc (nag_regsn_mult_linear_newyvar)
      - g02dkc (nag_regsn_mult_linear_tran_model)
      - g02dnc (nag_regsn_mult_linear_est_func)
      - g02eac (nag_all_regsn)
      - g02ecc (nag_cp_stat)
      - g02eec (nag_step_regsn)
      - g02fac (nag_regsn_std_resid_influence)
      - g02fcc (nag_durbin_watson_stat)
      - g02gac (nag_glm_normal)
      - g02gbc (nag_glm_binomial)
      - g02gcc (nag_glm_poisson)
      - g02gdc (nag_glm_gamma)
      - g02gkc (nag_glm_tran_model)
      - g02gnc (nag_glm_est_func)
      - g02hac (nag_robust_m_regsn_estim)
      - g02hbc (nag_robust_m_regsn_wts)
      - g02hdc (nag_robust_m_regsn_user_fn)
      - g02hfc (nag_robust_m_regsn_param_var)
      - g02hkc (nag_robust_corr_estim)
      - g02hlc (nag_robust_m_corr_user_fn)
      - g02hmc (nag_robust_m_corr_user_fn_no_derr)
    - g03 (Multivariate Methods)
    - g04 (Analysis of Variance)
    - g05 (Random Number Generators)
    - g07 (Univariate Estimation)
    - g08 (Nonparametric Statistics)
    - g10 (Smoothing in Statistics)
    - g11 (Contingency Table Analysis)
    - g12 (Survival Analysis)
    - g13 (Time Series Analysis)
    - h (Operations Research)
    - m01 (Sorting)
    - s (Approximations of Special Functions)
    - x01 (Mathematical Constants)
    - x02 (Machine Constants)
    - x04 (Input-Output Utilities)
  - Getting Started
  - Example Worksheet
エラーメッセージガイド
マニュアル
数学アプリ

NAG[g02hdc] NAG[nag_robust_m_regsn_user_fn] - Robust regression, compute regression with user-supplied functions and weights

Calling Sequence

g02hdc(chi, psi, psip0, beta, regtype, sigma_est, x, y, wgt, theta, k, sigma, rs, tol, eps, maxit, nitmon, nit, 'n'=n, 'm'=m, 'outfile'=outfile, 'comm'=comm, 'fail'=fail)

nag_robust_m_regsn_user_fn(. . .)

Parameters

chi - procedure;

If , chi must return the value of the weight function for a given value of its argument. The value of must be non-negative.

Its specification is:

chi(t, comm)

t - float;

On entry: the argument for which chi must be evaluated.

comm - table;

A Maple table, which should be generated using NAG[Nag_Comm], corresponding to the Nag_Comm structure.

The following fields may be accessed via NAG[GetOptions]:

p - Maple

Before calling nag_robust_m_regsn_user_fn (g02hdc) this field may be initialized for use by chi when called from nag_robust_m_regsn_user_fn (g02hdc).

chi is required only if , otherwise it can be specified as a pointer with 0 value.

psi - procedure;

psi must return the value of the weight function for a given value of its argument.

Its specification is:

psi(t, comm)

t - float;

On entry: the argument for which psi must be evaluated.

comm - table;

A Maple table, which should be generated using NAG[Nag_Comm], corresponding to the Nag_Comm structure.

The following fields may be accessed via NAG[GetOptions]:

p - Maple

Before calling nag_robust_m_regsn_user_fn (g02hdc) this field may be initialized for use by psi when called from nag_robust_m_regsn_user_fn (g02hdc).

psip0 - float;

On entry: the value of .

beta - float;

On entry: if , beta must specify the value of .

For Huber and Schweppe type regressions, is the 75th percentile of the standard Normal distribution (see g01fac (nag_deviates_normal)). For Mallows type regression is the solution to

where is the standard Normal cumulative distribution function (see s15abc (nag_cumul_normal)).

If , beta must specify the value of .

[[beta[2]=,(∫)chi(z)`ϕ`(z)dz , ],[,,],[,,],[beta[2]=,1/n(∑)w[i](∫)chi(z)`ϕ`(z)dz , ],[,,],[,,],[beta[2]=,1/n(∑)w[i]^2(∫)chi(z,/,w[i])`ϕ`(z)dz , ]]

where is the standard normal density, i.e., 1/(sqrt(2Pi))exp(-,1/2,x^2) .

If , beta is not referenced.

Constraint: if , . .

regtype - String;

On entry: determines the type of regression to be performed.

Huber type regression.

Mallows type regression.

Schweppe type regression.

sigma_est - String;

On entry: determines how is to be estimated.

is estimated by median absolute deviation of residuals.

is held constant at its initial value.

is estimated using the function.

x - Matrix(1..dim1, 1..m, datatype=float[8], order=order);

Note: this array may be supplied in Fortran_order or C_order , as specified by order. All array parameters must use a consistent order.

On entry: the values of the matrix, i.e., the independent variables. must contain the th element of , for ;

If , during calculations the elements of x will be transformed as described in Section [Description]. Before exit the inverse transformation will be applied. As a result there may be slight differences between the input x and the output x.

On exit: unchanged, except as described above.

y - Vector(1..n, datatype=float[8]);

On entry: the data values of the dependent variable.

must contain the value of for the th observation, for .

If , during calculations the elements of y will be transformed as described in Section [Description]. Before exit the inverse transformation will be applied. As a result there may be slight differences between the input y and the output y.

On exit: unchanged, except as described above.

wgt - Vector(1..n, datatype=float[8]);

On entry: the weight for the th observation, for .

If , during calculations elements of wgt will be transformed as described in Section [Description]. Before exit the inverse transformation will be applied. As a result there may be slight differences between the input wgt and the output wgt.

If , the th observation is not included in the analysis.

If , wgt is not referenced.

On exit: unchanged, except as described above.

theta - Vector(1..m, datatype=float[8]);

On entry: starting values of the argument vector . These may be obtained from least-squares regression. Alternatively if and or if and sigma approximately equals the standard deviation of the dependent variable, , then , for may provide reasonable starting values.

On exit: the M-estimate of , for .

k - assignable;

Note: On exit the variable k will have a value of type integer.

On exit: the column rank of the matrix .

sigma - assignable;

Note: On exit the variable sigma will have a value of type float.

On entry: a starting value for the estimation of . sigma should be approximately the standard deviation of the residuals from the model evaluated at the value of given by theta on entry.

Constraint: . .

On exit: the final estimate of if or the value assigned on entry if .

rs - Vector(1..n, datatype=float[8]);

On exit: the residuals from the model evaluated at final value of theta, i.e., rs contains the vector .

tol - float;

On entry: the relative precision for the final estimates. Convergence is assumed when both the relative change in the value of sigma and the relative change in the value of each element of theta are less than tol.

It is advisable for tol to be greater than .

Constraint: . .

eps - float;

On entry: a relative tolerance to be used to determine the rank of .

If or then machine precision will be used in place of tol.

A reasonable value for eps is where this value is possible.

maxit - integer;

On entry: the maximum number of iterations that should be used during the estimation.

A value of should be adequate for most uses.

Constraint: . .

nitmon - integer;

On entry: determines the amount of information that is printed on each iteration.

No information is printed.

On the first and every nitmon iterations the values of sigma, theta and the change in theta during the iteration are printed.

nit - assignable;

Note: On exit the variable nit will have a value of type integer.

On exit: the number of iterations that were used during the estimation.

'n'=n - integer; (optional)

On entry: , the number of observations.

Constraint: . .

'm'=m - integer; (optional)

On entry: , the number of independent variables.

Constraint: . .

'outfile'=outfile - character; (optional)

On entry: The name of a file to which intermediate or diagnostic output should be appended. If a value is not provided for this parameter then the behaviour of this routine is platform dependent. Usually all output will be suppressed, however on some platforms output will be produced and will be displayed in the Maple session.

'comm'=comm - table; (optional)

A Maple table, which should be generated using NAG[Nag_Comm], corresponding to the Nag_Comm structure.

'fail'=fail - table; (optional)

The NAG error argument, see the documentation for NagError.

Description

	Purpose
	nag_robust_m_regsn_user_fn (g02hdc) performs bounded influence regression (-estimates) using an iterative weighted least-squares algorithm.

Description

For the linear regression model

where is a vector of length of the dependent variable,

is a by matrix of independent variables of column rank ,

is a vector of length of unknown arguments,

and is a vector of length of unknown errors with var ,

nag_robust_m_regsn_user_fn (g02hdc) calculates the M-estimates given by the solution, , to the equation

(1)

where is the th residual, i.e., the th element of the vector ,

is a suitable weight function,

are suitable weights such as those that can be calculated by using output from g02hbc (nag_robust_m_regsn_wts),

and may be estimated at each iteration by the median absolute deviation of the residuals

or as the solution to

for a suitable weight function , where and are constants, chosen so that the estimator of is asymptotically unbiased if the errors, , have a Normal distribution. Alternatively may be held at a constant value.

The above describes the Schweppe type regression. If the are assumed to equal 1 for all , then Huber type regression is obtained. A third type, due to Mallows, replaces (1) by

This may be obtained by use of the transformations

[[w[i]^(*),:=sqrt(w[i])],[y[i]^(*),:=y[i]sqrt(w[i])],[x[ij]^(*),:=x[ij]sqrt(w[i]) j=1,2,..,m]]

(see Marazzi (1987b)).

The calculation of the estimates of can be formulated as an iterative weighted least-squares problem with a diagonal weight matrix given by

G[ii]={[[(psi(r[i],/,(sigma,w[i])))/((r[i],/,(sigma,w[i]))) ,r[i]<>0],[,],[psi'(0) ,r[i]=0 ]]

The value of at each iteration is given by the weighted least-squares regression of on . This is carried out by first transforming the and by

[[(y˜)[i],=y[i]sqrt(G[ii])],[(x˜)[ij],=x[ij]sqrt(G[ii]) j=1,2,..,m]]

and then using a least squares solver. If is of full column rank then an orthogonal-triangular (QR) decomposition is used; if not, a singular value decomposition is used.

Observations with zero or negative weights are not included in the solution.

Note: there is no explicit provision in the function for a constant term in the regression model. However, the addition of a dummy variable whose value is for all observations will produce a value of corresponding to the usual constant term.

nag_robust_m_regsn_user_fn (g02hdc) is based on routines in ROBETH, see Marazzi (1987b).

Error Indicators and Warnings

"NE_ALLOC_FAIL"

Dynamic memory allocation failed.

"NE_BAD_PARAM"

On entry, argument had an illegal value.

"NE_CHI"

Value given by chi function : .

"NE_CONVERGENCE_SOL"

Iterations to solve weighted least squares equations failed to converge.

"NE_CONVERGENCE_THETA"

Iterations to calculate estimates of theta failed to converge in maxit iterations: .

"NE_FULL_RANK"

Weighted least squares equations not of full rank: rank .

"NE_INT"

On entry, . Constraint: .

"NE_INT_2"

On entry, , . Constraint: .

On entry, : , .

"NE_INTERNAL_ERROR"

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please consult NAG for assistance.

"NE_NOT_CLOSE_FILE"

Cannot close file .

"NE_NOT_WRITE_FILE"

Cannot open file for writing.

"NE_REAL"

On entry, . Constraint: .

"NE_ZERO_DF"

Value of : , .

"NE_ZERO_VALUE"

Estimated value of sigma is zero.

	Accuracy
	The accuracy of the results is controlled by tol.

Further Comments

In cases when it is important for the value of sigma to be of a reasonable magnitude. Too small a value may cause too many of the winsorized residuals, i.e., , to be zero, which will lead to convergence problems and may trigger the the error "NE_FULL_RANK" is raised error.

By suitable choice of the functions chi and psi this function may be used for other applications of iterative weighted least-squares.

For the variance-covariance matrix of see g02hfc (nag_robust_m_regsn_param_var).

Examples

chi := proc(t, comm)
  local ps:
  if (abs(t) < 1.5) then
    ps := t:
  else
    ps := 1.5:
  end if:
  return (ps^2/2.0):
end proc:
psi := proc(t, comm)
  if (t <= -1.5) then
    -1.5:
  elif (t >= 1.5) then
    1.5:
  else
    t:
  end if:
end proc:
psip0 := 1:
beta := 0.1443849979905463:
regtype := "Nag_SchweppeReg":
sigma_est := "Nag_SigmaChi":
n := 5:
m := 3:
sigma := 1:
tol := 5e-05:
eps := 5e-06:
maxit := 50:
nitmon := 0:
comm := NAG:-Nag_Comm():
x := Matrix([[1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1], [1, 0, 3]], datatype=float[8]):
y := Vector([10.5, 11.3, 12.6, 13.4, 17.1], datatype=float[8]):
wgt := Vector([0.4039, 0.5012, 0.4039, 0.5012, 0.3862], datatype=float[8]):
theta := Vector([0, 0, 0], datatype=float[8]):
rs := Vector(5, datatype=float[8]):
NAG:-g02hdc(chi, psi, psip0, beta, regtype, sigma_est, x, y, wgt, theta, k, sigma, rs, tol, eps, maxit, nitmon, nit, 'n' = n, 'm' = m, 'comm' = comm):