% % Construct multivariate polynomial response surface for exponential model % % Generate 300 training points %copyfile('PolyfitnTools.mltbx') % matlab.addons.toolbox.installToolbox('PolyfitnTools.mltbx') % %% copyfile('PolyfitnTools.mltbx'); %unzip('PolyfitnTools.mltbx'); % % Compute the solution at n_train = 300 training points and add sigma=0.1 random noise. m = 2; n_train = 300; q_train = -1+2*rand(n_train,m); f_train = exp(0.7*q_train(:,1)+0.3*q_train(:,2))+normrnd(0,0.01,n_train,1); % % Compute multiviariate polynomials of orders 1-5 and compute AIC and BIC. % sigma=.1; for i = 1:5 p = polyfitn(q_train,f_train,i); k = length(p.Coefficients); ssq = sum((polyvaln(p,q_train)-f_train).^2); % Sum of squares L = (1/((2*pi*sigma^2)^(n_train/2)))*exp(-ssq/(2*sigma^2)); % Gaussian likelihood aic(i) = 2*k-2*log(L); bic(i) = -2*log(L)+k*log(n_train); end order = find(aic == min(aic)); % Find minimium AIC value % % Construct response surface with appropriate polynomial % p = polyfitn(q_train,f_train,order(1)); qplot1 = linspace(-1,1,100); qplot2 = linspace(-1,1,100); [z1 z2] = meshgrid(qplot1,qplot2); pairs = [z1(:) z2(:)]; resp = reshape(polyvaln(p,pairs),100,100); % % Generate 300 testing points % n_test= 300; q_test= -1+2*rand(n_test,m); f_test = exp(0.7*q_test(:,1)+0.3*q_test(:,2))+normrnd(0,sigma,n_test,1); % % Compute RMSE % resp_test = polyvaln(p,q_test); rmse = sqrt((1/n_test)*sum((f_test-resp_test).^2)); % % Plot response surface with testing data % figure(3) plot3(q_test(:,1),q_test(:,2),f_test,'ob','Linewidth',2) hold on surf(qplot1, qplot2, resp); hold off