%% example of nlinfit using file name = inv_soln.m
%This program can be used as a base for most nonlinear regression OLS
%using nlinfit
%% Housekeeping
%
clear %all; % Clear the workspace.
close all; % Close all figures.
format compact
%% Read in data
data =xlsread('exp_data.xls');
x=data(:,1);
yobs=data(:,2);
%% Initial parameter guesses
yo=100;
k=0.6;
beta0(1)=yo; %initial guess yo
beta0(2)=k; %initial guess kr
%% X' = scaled sensitivity coefficients using forward-difference
% This is a forward problem with known approximate parameters
%Xp=SSC(beta0,x,@forderdiff);
Xp=SSC(beta0,x,@forderexp);
title('scaled sensitivity coefficients using initial guesses')
type SSC.m
type forderdiff.m
%% nlinfit returns parameters, residuals, Jacobian (sensitivity %coefficient matrix),
%covariance matrix, and mean square error. ode45 is solved many times iteratively
[beta,resids,J,COVB,mse] = nlinfit(x,yobs,@forderdiff,beta0);
%[beta,resids,J,COVB,mse] = nlinfit(x,yobs,@forderexp,beta0);
beta
rmse=sqrt(mse)
condX=cond(J)
detXTX=det(J'*J)
% type forderdiff.m %uncomment if you wish to print the code of the function
%% confidence intervals for parameters
ci=nlparci(beta, resids,J)
%% R is the correlation matrix for the parameters, sigma is the standard error vector
[R,sigma]=corrcov(COVB);
R
sigma
relerr=sigma./beta' %relative error (coefficient of variance) for each parameter
%% Confidence and prediction intervals for the dependent variable
%nonlinear regression confidence intervals-- 'on' means simultaneous
%bounds; 'off' is for nonsimultaneous bounds; must use 'curve' for
%regression line, 'observation' for prediction interval
% [ypred, delta] = nlpredci('forderdiff',x,beta,resids,J,0.05,'on','curve'); %confidence band for regression line
% [ypred, deltaob] =nlpredci('forderdiff',x,beta,resids,J,0.05,'on','observation');%prediction band for individual points
[ypred, delta] = nlpredci('forderexp',x,beta,resids,J,0.05,'on','curve'); %confidence band for regression line
[ypred, deltaob] =nlpredci('forderexp',x,beta,resids,J,0.05,'on','observation');%prediction band for individual points
yspan=range(ypred)% total span of ypred
relrmse=rmse/yspan % ratio of rmse vs. yspan
%simultaneous confidence bands for regression line
CBu=ypred+delta;
CBl=ypred-delta;
%simultaneous prediction bands for regression line
PBu=ypred+deltaob;
PBl=ypred-deltaob;
%% Output--ypred and yobs vs. t
figure
hold on
h1(1)=plot(x,ypred,'-','linewidth',3); %predicted y values
h1(2)=plot(x,yobs,'square', 'Markerfacecolor', 'r');
xlabel('time (min)','fontsize',16,'fontweight','bold')
ylabel('y','fontsize',16,'fontweight','bold')
grid on
%% Output --CIs and PIs
%plot Cobs, Cpred line, confidence band for regression line
h1(3) = plot(x,CBu,'--g','LineWidth',2);
plot(x,CBl,'--g','LineWidth',2);
%plot prediction band for regression line
h1(4) = plot(x,PBu,'-.','LineWidth',2);
plot(x,PBl,'-.','LineWidth',2);
legend(h1,'ypred','yobs','CB','PB')
%% residual scatter plot
figure
hold on
plot(x, resids, 'square','Markerfacecolor', 'b');
YLine = [0 0];
XLine = [0 max(x)];
plot (XLine, YLine,'R'); %plot a straight red line at zero
ylabel('Observed y - Predicted y','fontsize',16,'fontweight','bold')
xlabel('time (min)','fontsize',16,'fontweight','bold')
grid on
%% residuals histogram--same as dfittool, but no curve fit here
[n1, xout] = hist(resids,10); %10 is the number of bins
figure
hold on
set(gca, 'fontsize',14,'fontweight','bold');
bar(xout, n1) % plots the histogram
xlabel('Y_{observed} - Y_{predicted}','fontsize',16,'fontweight','bold')
ylabel('Frequency','fontsize',16,'fontweight','bold')
%% scaled sensitivity coefficients using final estimated parameters
% This is a double-check to make sure X' has not changed much
Xp=SSC(beta,x,@forderexp);
%Xp=SSC(beta,x,@forderexp);
title('scaled sensitivity coefficients using beta estimates')
nlinfit :Nonlinear least-squares regression.
BETA = nlinfit(X,Y,MODELFUN,BETA0) estimates the coefficients of a
nonlinear regression function, using least squares estimation. Y is a
vector of response (dependent variable) values. Typically, X is a
design matrix of predictor (independent variable) values, with one row
for each value in Y and one column for each coefficient. However, X
may be any array that MODELFUN is prepared to accept. MODELFUN is a
function, specified using @, that accepts two arguments, a coefficient
vector and the array X, and returns a vector of fitted Y values. BETA0is a vector containing initial values for the coefficients.
[BETA,R,J,COVB,MSE] = nlinfit(X,Y,MODELFUN,BETA0) returns the fitted
coefficients BETA, the residuals R, the Jacobian J of MODELFUN, the
estimated covariance matrix COVB for the fitted coefficients, and an
estimate MSE of the variance of the error term. You can use these
outputs with NLPREDCI to produce confidence intervals for predictions,
and with NLPARCI to produce confidence intervals for the estimated
coefficients. If you use a robust option (see below), you must use
COVB and may need MSE as input to NLPREDCI or NLPARCI to insure that
the confidence intervals take the robust fit properly into account.
转载本文请联系原作者获取授权,同时请注明本文来自陆宇振科学网博客。
链接地址:http://blog.sciencenet.cn/blog-578676-836767.html
上一篇:Solving Second-oder Differential Equation with Matlab ode45
下一篇:Computing Jocobian and Scaled Sensitivity Coefficients