Regression Model Graph(plot), Deflection Problem(sharp point or cusp)

I implemented the first, second, and third regression models for the data and showed them using the plot function.

But, for the 3rd regression model, the graph was bent, such as the part in red.

How can we solve this problem?

[KSSV]

Why do you think it is a problem? When you increase the order, polynomial is increased and it will get close to the given points.

[Sanghyuk Kim]

yeah i know.

I think that the graph of the 3rd order regression model should be a smooth curve.

However, as you can see, the graph is bent. I think this is the problem.

[Scott MacKenzie]

It might help if you post the data and the code that generated the plot you are seeing

[Sanghyuk Kim]

load accidents

x = hwydata(:,6); %Population of states

y = hwydata(:,4); %Accidents per state

scatter(x,y,'filled')

X = [ones(size(x)) x];

b = regress(y,X)

hold on

YFIT = b(1) + b(2)*x;

plot(x,YFIT)

X2 = [ones(size(x)) x x.^2];

b2 = regress(y,X2);

YFIT2 = b2(1) + b2(2)*x + b2(3)*x.^2 ;

plot(x,YFIT2)

X3 = [ones(size(x)) x x.^2 x.^3];

b3 = regress(y,X3);

YFIT3 = b3(1) + b3(2)*x + b3(3)*x.^2 + b3(4)*x.^3 ;

plot(x,YFIT3)

xlabel('Registered vehicles[thousands](x)');

ylabel('Accidents per state(y)');

plot1_legend=legend('Data','1st Order','2nd Order', '3rd Order')

hold off

The data is Accidents data set!

[Sanghyuk Kim]

I think a 3rd regression model is like a 3rd polynomial, so there should be no bends.

[Sanghyuk Kim]

3rd regression model is

YFIT3 = 201.6374 + 0.0815*x + (9.7323e-6)*x.^2 - (2.6249e-10)*x.^3

My guess is that the error is probably due to the coefficient being too small.

how to fix it?

Answer 1

[Scott MacKenzie]

You are getting the sharp bend because you are calculating YFIT using the x sample points (which are sparce at the inflexion point). To see the curvature in the model, you need to calculate YFIT using a set of x query points (xq) that are finely spaced over the full range of interest (including about the inflexion point):

load accidents
x = hwydata(:,6); % population of states
y = hwydata(:,4); % accidents per state
scatter(x,y,'filled')
hold on;

% use xq for plotting (to expose curvature in models)  
xq = linspace(min(x), max(x));

X = [ones(size(x)) x];
b = regress(y,X)   
hold on
YFIT = b(1) + b(2)*xq; % use xq, not x
plot(xq,YFIT);
hold on;

X2 = [ones(size(x)) x x.^2];
b2 = regress(y,X2);
YFIT2 = b2(1) + b2(2)*xq + b2(3)*xq.^2 ; % use xq, not x
plot(xq,YFIT2);

X3 = [ones(size(x)) x x.^2 x.^3];
b3 = regress(y,X3);
YFIT3 =  b3(1) + b3(2)*xq + b3(3)*xq.^2 + b3(4)*xq.^3 ; % use xq, not x
plot(xq,YFIT3) ;

xlabel('Registered vehicles[thousands](x)');
ylabel('Accidents per state(y)');
legend('Data','1st Order','2nd Order', '3rd Order', 'location', 'southeast');

BTW, consider using polyfit (along with polyval) instead of regress. It's a bit easier overall:

pf1 = polyfit(x, y, 1)
yfit1 = polyval(pf1, xq)
plot (xq, yfit1)

pf2 = polyfit(x, y, 2)
yfit2 = polyval(pf2, xq)
plot (xq, yfit2)

pf3 = polyfit(x, y, 3)
yfit3 = polyval(pf3, xq)
plot (xq, yfit3)