Main Content

nlpredci

非線形回帰予測の信頼区間

説明

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Covar',CovB) は、入力値 X について、非線形回帰モデル modelfun に対する予測値 Ypred と 95% の信頼区間の半値幅 delta を返します。nlpredci を呼び出す前に modelfun を近似するため、nlinfit を使用して、推定された係数値 beta、残差 R および分散共分散行列 CovB を取得してください。

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Covar',CovB,Name,Value) は、1 つ以上の名前と値のペアの引数によって指定された追加オプションを使用します。

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Jacobian',J) は、入力値 X について、非線形回帰モデル modelfun に対する予測値 Ypred と 95% の信頼区間の半値幅 delta を返します。nlpredci を呼び出す前に、modelfun を近似するために nlinfit を使用し、推定された係数値 beta、残差 R およびヤコビ行列 J を取得してください。

nlinfit でロバストなオプションを使用する場合、Jacobian 構文の代わりに Covar 構文を使用してください。分散共分散行列 CovB は、ロバスト近似を正しく考慮しなければなりません。

[Ypred,delta] = nlpredci(modelfun,X,beta,R,'Jacobian',J,Name,Value) は、1 つ以上の名前と値のペアの引数によって指定された追加オプションを使用します。

すべて折りたたむ

標本データを読み込みます。

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

beta0 にある初期値を使用して、レート データに Hougen-Watson モデルを当てはめます。

[beta,R,J] = nlinfit(X,y,@hougen,beta0);

平均反応物レベルにおける、曲線の値に対する予測応答および 95% の信頼区間の半値幅を取得します。

[ypred,delta] = nlpredci(@hougen,mean(X),beta,R,'Jacobian',J)
ypred = 5.4622
delta = 0.1921

曲線の値について 95% の信頼区間を計算します。

[ypred-delta,ypred+delta]
ans = 1×2

    5.2702    5.6543

標本データを読み込みます。

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

beta0 にある初期値を使用して、レート データに Hougen-Watson モデルを当てはめます。

[beta,R,J] = nlinfit(X,y,@hougen,beta0);

反応物レベル [100,100,100] における、新しい観測の予測応答および 95% の予測区間の半値幅を取得します。

[ypred,delta] = nlpredci(@hougen,[100,100,100],beta,R,'Jacobian',J,...
                         'PredOpt','observation')
ypred = 1.8346
delta = 0.5101

新しい観測の 95% の予測区間を計算します。

[ypred-delta,ypred+delta]
ans = 1×2

    1.3245    2.3447

非線形回帰モデル y=b1+b2exp{b3x}+ϵ から標本データを生成します。ここで、b1b2 および b3 は係数です。誤差項は、平均が 0、標準偏差が 0.5 の正規分布に従います。

modelfun = @(b,x)(b(1)+b(2)*exp(-b(3)*x));

rng('default') % for reproducibility
b = [1;3;2];
x = exprnd(2,100,1);
y = modelfun(b,x) + normrnd(0,0.5,100,1);

ロバスト近似オプションを使用して非線形モデルを当てはめます。

opts = statset('nlinfit');
opts.RobustWgtFun = 'bisquare';
beta0 = [2;2;2];
[beta,R,J,CovB,MSE] = nlinfit(x,y,modelfun,beta0,opts);

近似した回帰モデルと 95% の同時信頼限界をプロットします。

xrange = min(x):.01:max(x);
[ypred,delta] = nlpredci(modelfun,xrange,beta,R,'Covar',CovB,...
                         'MSE',MSE,'SimOpt','on');
lower = ypred - delta;
upper = ypred + delta;

figure()
plot(x,y,'ko') % observed data
hold on
plot(xrange,ypred,'k','LineWidth',2)
plot(xrange,[lower;upper],'r--','LineWidth',1.5)

Figure contains an axes object. The axes object contains 4 objects of type line. One or more of the lines displays its values using only markers

標本データを読み込みます。

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

観測の重みの関数ハンドルを指定した後、指定した観測の重み関数を使用してレート データに Hougen-Watson モデルを当てはめます。

a = 1; b = 1;
weights = @(yhat) 1./((a + b*abs(yhat)).^2);
[beta,R,J,CovB] = nlinfit(X,y,@hougen,beta0,'Weights',weights);

観測の重み関数を使用して、反応物レベル [100,100,100] における新しい観測の 95% の信頼区間を計算します。

[ypred,delta] = nlpredci(@hougen,[100,100,100],beta,R,'Jacobian',J,...
                         'PredOpt','observation','Weights',weights);
[ypred-delta,ypred+delta]
ans = 1×2

    1.5264    2.1033

標本データを読み込みます。

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

組み合わせた誤差分散モデルを使用して、レート データに Hougen-Watson モデルを当てはめます。

[beta,R,J,CovB,MSE,S] = nlinfit(X,y,@hougen,beta0,'ErrorModel','combined');

近似した誤差分散モデルを使用して、反応物レベル [100,100,100] における新しい観測値の 95% の予測区間を計算します。

[ypred,delta] = nlpredci(@hougen,[100,100,100],beta,R,'Jacobian',J,...
                         'PredOpt','observation','ErrorModelInfo',S);
[ypred-delta,ypred+delta]
ans = 1×2

    1.3245    2.3447

入力引数

すべて折りたたむ

関数ハンドルとして指定される、非線形回帰モデル関数。modelfun は、係数ベクトルと配列 X の 2 つの入力引数を順序どおりに受け入れ、近似した応答値のベクトルを返さなければなりません。

たとえば、非線形回帰関数 hougen を指定するには、関数ハンドル @hougen を使用します。

データ型: function_handle

行列として指定される、予測の入力値。nlpredci は、X の各行にある共変量を予測します。X にはモデルの各係数に対して 1 つの列が必要です。

データ型: single | double

推定回帰係数。nlinfit の以前の呼び出しで返された近似係数のベクトルとして指定します。

データ型: single | double

近似された modelfun の残差。nlinfit の以前の呼び出しで返された残差のベクトルとして指定します。

近似された係数 beta の推定分散共分散行列。nlinfit の以前の呼び出しで返された分散共分散行列として指定します。

非線形回帰モデル modelfun の推定されたヤコビアン。nlinfit の以前の呼び出しで返されたヤコビ行列として指定します。

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後ろにする必要がありますが、ペアの順序は関係ありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: 'Alpha',0.1,'PredOpt','observation' — 新しい観測の 90% の予測区間

信頼区間の有意水準。'Alpha' と、(0,1) の範囲内のスカラー値で構成されるコンマ区切りのペアとして指定します。Alpha の値が α の場合、nlpredci は 100×(1–α)% の信頼水準の区間を返します。

既定の信頼水準は 95% (α = 0.05) です。

例: 'Alpha',0.1

データ型: single | double

誤差モデルの近似に関する情報。'ErrorModelInfo' と、nlinfit の以前の呼び出しで返された構造体で構成されるコンマ区切りのペアとして指定します。

ErrorModelInfo は、PredOpt の値が 'observation' の場合にのみ、返される予測区間に影響します。ErrorModelInfo を使用しない場合、nlpredci は誤差分散モデルが 'constant' であると仮定します。

nlinfit によって返される誤差モデル構造体には、以下のフィールドがあります。

ErrorModel選択された誤差モデル
ErrorParameters推定された誤差パラメーター
ErrorVarianceN 行 p 列の行列 X を受け入れ、推定された誤差モデルを使用して誤差分散の N 行 1 列のベクトルを返す関数ハンドル。
MSE平均二乗誤差
ScheffeSimPred推定された誤差モデルを使用する場合の、同時予測区間の Scheffé パラメーター
WeightFunctionnlinfit で以前にカスタムの重み関数を使用している場合は値が true の論理値
FixedWeightsnlinfit で以前に固定された重みを使用している場合は値が true の論理値
RobustWeightFunctionnlinfit で以前にロバスト近似を使用した場合は値が true の論理値

近似された非線形回帰モデルの平均二乗誤差 (MSE)。'MSE' と、nlinfit の以前の呼び出しで返された MSE 値で構成されるコンマ区切りのペアとして指定します。

nlinfit でロバスト オプションを使用する場合、ロバスト近似が正しく考慮されるよう、新しい観測の予測時には MSE を指定しなければなりません。MSE を指定しない場合、nlpredci は残差 R から MSE を計算し、ロバスト近似を考慮に入れません。

たとえば、nlinfit から返された MSE の値が mse である場合、'MSE',mse と指定できます。

データ型: single | double

計算対象の予測区間。'PredOpt''curve' または 'observation' で構成されるコンマ区切りのペアとして指定します。

  • 値を 'curve' に指定すると、nlpredci は観測値 X において推定された曲線 (関数値) の信頼区間を返します。

  • 値を 'observation' に指定すると、nlpredciX における新しい観測値の信頼区間を返します。

nlinfit でロバストなオプションを使用した後で 'observation' を指定する場合、平均二乗誤差のロバスト推定を指定するには、MSE の値も指定しなければなりません。

例: 'PredOpt','observation'

同時区間を指定するためのインジケーター。'SimOpt' と、'off' または 'on' で構成される、コンマ区切りペアとして指定します。非同時区間を計算する場合は 'off'、同時区間を計算する場合は 'on' の値を使用します。

観測値の重み。'Weights' と、正のスカラー値のベクトルまたは関数ハンドルで構成されるコンマ区切りペアとして指定します。既定の設定は重みなしです。

  • 重みのベクトルを指定する場合、ベクトルには X にある観測値 (行) と同じ数の要素を含めなければなりません。

  • 重みの関数のハンドルとして指定する場合、その関数ハンドルは予測される応答値のベクトルを入力として受け入れ、出力として正の実数重みのベクトルを返さなければなりません。

重み W に対して nlpredci は観測 i における誤差分散を mse*(1/W(i)) により推定します。mse は、MSE を使用して求めた平均二乗誤差値です。

例: 'Weights',@WFun

データ型: double | single | function_handle

出力引数

すべて折りたたむ

予測された応答であり、X と同じ行数のベクトルとして返されます。

X と同じ行数のベクトルとして返される、信頼区間の半値幅。既定での delta には、X の観測値における modelfun に関する 95% の非同時信頼区間が含まれます。信頼区間の下限と上限はそれぞれ Ypred-delta および Ypred+delta として計算できます。

'PredOpt' の値が 'observation' である場合、delta には X の値における新しい観測値の予測区間の半値幅が含まれます。

詳細

すべて折りたたむ

推定可能な予測値の信頼区間

推定されたモデルのヤコビ行列がフル ランクでない場合、すべての予測点において妥当な信頼区間を作成できない可能性があります。その場合でも nlpredci はすべての "推定可能" な予測点について信頼区間を作成しようと試みます。

たとえば、次の計画行列内の点で線形関数 f(xi,β)=β1xi1+β2xi2+β3xi3 を当てはめるとします。

X=(110110110101101101).

X の値で推定されるヤコビ行列は計画行列自体になるので J=X. です。したがって、このヤコビ行列はフル ランクではありません。

rng('default') % For reproducibility
y = randn(6,1);

linfun = @(b,x) x*b;
beta0 = [1;1;1];
X = [repmat([1 1 0],3,1); repmat([1 0 1],3,1)];   

[beta,R,J]  = nlinfit(X,y,linfun,beta0);
Warning: The Jacobian at the solution is ill-conditioned, and
some model parameters may not be estimated well (they are not
identifiable).  Use caution in making predictions. 
> In nlinfit at 283 

この例では、nlpredci は次の線形関係を満たす点における予測区間のみを計算できます。

xi1=xi2+xi3.

識別できない点で予測の信頼区間を計算しようとすると、nlpredci は対応する信頼区間の半値幅について NaN を返します。

xpred = [1 1 1;0 1 -1;2 1 1];
[ypred,delta] = nlpredci(linfun,xpred,beta,R,'Jacobian',J)
ypred =

   -0.0035
    0.0798
   -0.0047


delta =

       NaN
    3.8102
    3.8102
delta の最初の要素が NaN になるのは、xpred の最初の行が必要な線形依存関係を満たさず、推定可能な対比ではないためです。

ヒント

  • 複雑なパラメーターまたはデータの信頼区間を計算するには、問題をその実数部と虚数部に分割する必要があります。nlinfit を呼び出すときは、次の手順を実行します。

    1. パラメーター ベクトル beta を、元のパラメーター ベクトルの実数部と虚数部の連結として定義します。

    2. 応答ベクトル Y の実数部と虚数部を単一のベクトルとして連結します。

    3. モデル関数 modelfun を変更し、X と純粋な実数パラメーター ベクトルを受け入れて、近似された値の実数部と虚数部の連結が返されるようにします。

    このように問題を定式化して、nlinfit は実際の推定値を計算し、信頼区間は実行可能になります。

アルゴリズム

  • nlpredci は、残差 R またはヤコビ行列 JNaN 値を欠損値として認識し、対応する観測を無視します。

  • ヤコビ行列 J の列がフル ランクでない場合、モデル パラメーターの一部が認識されないことがあります。その場合、nlpredci推定された予測に対する信頼区間の作成を試みますが、パラメーターが認識されない場合は NaN を返します。

参照

[1] Lane, T. P. and W. H. DuMouchel. “Simultaneous Confidence Intervals in Multiple Regression.” The American Statistician. Vol. 48, No. 4, 1994, pp. 315–321.

[2] Seber, G. A. F., and C. J. Wild. Nonlinear Regression. Hoboken, NJ: Wiley-Interscience, 2003.

バージョン履歴

R2006a より前に導入