# Fit a smooth approximation line

I have a simulated data set and I want to fit a "smooth approximation line" like the image I have provided. 1- Is it possible to do this in Excel or Matlab? May I have a pointer on how to do it?

2- How we find how much of the data is within a distance of 0.15 of the line? Or do we basically fit a line in a way that 90% of the points fall in 0.15 range?

• I have updated my answer a third time. Without being sure what the model is, we are left to guess at how to emulate it. – Glen_b -Reinstate Monica Mar 29 '13 at 5:49
• I have finally addressed the two parts of your question directly; see part 6 of this opus. More broadly, you should probably consult a good regression text. – Glen_b -Reinstate Monica Mar 30 '13 at 0:07
• I added some more discussion of how you do the interval stuff ... and we see maybe where the 0.15 might have come from at the end. – Glen_b -Reinstate Monica Mar 30 '13 at 7:17

There are many ways to smooth. Depending on how much effort you're prepared to go to, any of can be done in matlab (and probably already exist somewhere) and Excel (but you might have to do them yourself).

[If you were asking about R I'd show you three or four ways you might smooth data.]

That "smooth approximation" in the image just looks quadratic to me, but it might be say a natural spline or something similar.

Edit:

On reading around it looks like it's possible that there's some particular functional form, but I have had no success finding out what it is. It's not possible to tell what they did from the information I can see.

If you know what functional form was actually fitted, please state it in your question.

If you don't know, perhaps the best thing to do is write to the paper's authors.

Details seem to be quite hard to find.

Edit 2:

Well, I've found the data, which is a start.

The data are here, specifically the mixed work traffic subset of that data; there are two variables, v and rho, representing density and velocity. The plot you have is (as I eventually discovered from reading many papers) the product of v and rho ("flow") plotted against rho (density).

Here's what the originally speculated "quadratic" looks like (in this case, constrained to go through the origin: As you see, it's very similar to the plot you have, but not quite identical. The difference seems to be that I simply used least squares regression, whereas the paper your plot is from:

Davidich, M. and G. Köster (2012),
Towards automatic and robust adjustment of human behavioral parameters in a pedestrian stream model to measured data
Safety Science, Vol., No. 5, June, p. 1253–1260

seems to suggest that perhaps they used a different approach to obtain their fit (though apart from at the very end the fit seems to be essentially identical). The difference may even just be rounding.

Nevertheless I'm satisfied their plot is very likely actually quadratic in this case.

some R code for anyone interested in playing with it:

rho <- c(0.68, 0.15, 0.2, 0.3, 0.35, 0.4, 0.48, 0.55, 0.6, 0.9, 1.1,
1.4, 1.5, 1.6, 1.65, 1.7, 2, 2.1, 2.5, 2.6, 1.3, 1.3, 1.4, 1.5,
1.55, 1.64, 1.7, 1.9, 0.55, 0.95, 1.2, 1.2, 1.4, 1.5, 1.55, 1.85,
0.65, 0.8, 0.8, 0.4, 0.5, 0.4, 0.58, 0.4, 0.75, 0.45, 0.9, 0.65, 0.75)

v <- c(1.34, 1.65, 1.45, 1.5, 1.45, 1.4, 1.4, 1.35, 1.3, 1.1, 0.95,
1.1, 1.1, 0.9, 0.85, 0.8, 0.77, 0.6, 0.7, 0.6, 1.1, 1.2, 1, 1.25,
1.15, 1.06, 1, 0.85, 1.45, 1.25, 1.15, 1.1, 0.9, 1, 0.95, 0.85,
1.3, 1.3, 1.4, 1.35, 1.4, 1.45, 1.3, 1.5, 1.28, 1.55, 1.2, 1.35, 1.5)

J <- rho*v
plot(J~rho,xlim=c(0,2.6),ylim=c(0,2),col=4,cex=.7,pch=16)
mwt.lm <- lm(J~rho+I(rho^2)-1)

xx <- seq(0,2.6,0.05)
pr <- predict.lm(mwt.lm,newdata=data.frame(rho=xx))
lines(pr~xx)


Edit 3:

Well, they're not quite as close as I thought - here they are superimposed (you can only see my points because my dots are slightly bigger). The thick curve is the original, the thinner black curve is my least-squares-through-the-origin-quadratic above: The fit from the paper might still be quadratic, but if it is the peak is noticeably further to the right than mine (mine is about 2.055 while theirs is somewhere around 2.2). I suspect they might have done something else though.

Edit 4:

I notice that the velocity-density diagram looks vaguely negative exponential in form, which suggests the possibility of a Hoerl-type curve ($$E(J|\rho) = a\rho\exp{(-b\rho)}$$). The fit is in green (the red dashed one is my earlier quadratic). Their fitted curve seems to be in between the two.

I tend to lean back toward my suggestion of some form of spline, probably a natural or semi-natural cubic spline, but there may well be some specific functional form involved

One might speculate at length; I suspect your best bet is to ask the authors what they actually did.

Here's one way to fit a curve like that, done in R:

hoerl.lm <- lm(log(J)~rho+log(rho))
prh <- predict.lm(hoerl.lm,newdata=data.frame(rho=xx))
lines(exp(prh)~xx,col=3)


(assumes the previous data and plot code has been run)

Edit 5:

General advice on fitting models to this sort of data:

From the data sets I've seen $$v$$ vs $$\rho$$ tends to be very nearly linear over a wide range of values, and it tends to have more nearly constant variance.

As a result, I strongly suggest modelling $$v$$ as a function of $$\rho$$ and then if you want a model for $$J = \rho v$$, you multiply the model for $$v$$ by $$\rho$$ (taking care that if the original model for $$v$$ had constant variance, the scale of the random variation about the model for $$J$$ will now be proportional to $$\rho$$.

If you fit a linear model on the $$v$$ scale, you'll get a quadratic-through-the-origin for $$J$$.

So for example, for this data, this is what $$v$$ vs $$\rho$$ looks like: Clearly quite close to linear. I tried a variety of smooths - natural splines with various degrees of freedom and loess type smooths with various degrees and spans, but they are all pretty close to linear, so when converted to the $$J$$ vs $$\rho$$ scale they mostly look fairly close to the quadratic fit above. Some were flatter on the right side, but none came any closer to reproducing the original fit.

Edit 6:

A quadratic fit which better reproduces the characteristics of the fit in the paper would be:

$$J = 1.525 \rho -0.356 \rho^2\,$$ .

I don't think the fit is actually quadratic, but quadratics somewhere in this region of parameter space do look a lot like the curve in the paper; it's pretty close.

Now to the outstanding original questions

I'll describe some details for some specific basic cases.

1) how to fit a least squares quadratic directly to the $$J$$ vs $$\rho$$ data.

a) You need to fit rho and rho*rho as predictors. In Excel this is best done with the Data Analysis toolpack; in Matlab the regression command is of the form [b,bint,r,rint,stats] = regress(y,X)

• but to have the quadratic you would put $$\rho$$ and $$\rho^2$$ in as columns of X and omit the intercept by leaving out the constant column of 1's.

b) if you want it to go through the origin (and you definitely do), then you need to omit the intercept term.

How we find how much of the data is in the 0.15 range of the line? Or we basically fit a line in a way that 90% of the pints fall in 0.15 range?

There are two possibilities:

a) that 0.15 is of particular interest, and you want to know how much of the data is within that distance of the line. You can estimate that in a nonparametric way simply by counting to get a proportion of points within that distance. If the variance is constant, would also be an estimate of the proportion in the population and hence in another sample - otherwise you're stuck with worrying about the relative proportions at each density in the samples or in sample and population. You can do it parametrically using normal assumptions as well, by finding the proportion of a $$t_{n-p}$$ distribution between $$-0.15/s$$ and $$0.15/s$$ where $$s$$ is the regression's standard deviation of the residuals (with a divisor of $$n-p$$). Note that $$p$$ is the number of predictors in the model (including an intercfept if you have one).

* (95th percentile for a 90% of the data, because you'll lose 5% out each end)

or - more likely -

b) that 90% is of particular interest and you want to see what distance from the line you need to go to get 90%. You can do it nonparametrically by counting out to an interval containing 90% of the data (easiest done by sorting the absolute residuals and finding their 90th percentile), or you can do it parametrically under a normal assumption by multiplying the regression's $$s$$ by the 95th percentile* of a $$t_{n-p}$$ distribution (for your sample size it's 1.678).

However, these "90% within 0.15" calculations - if you want to get much inferential value out of them at least - rely on the assumption of constant spread. You simply don't have that; the relationship between $$v$$ and $$\rho$$ seems to have nearly constant spread here, and so the spread of $$J$$ values around its curve is proportional to $$\rho$$. It makes more sense to describe your intervals on that scale and then multiply by $$\rho$$ to produce intervals on the $$J$$-scale.

--

Now to model it on the $$v$$ vs $$\rho$$ scale, you regress $$v$$ on $$\rho$$ with the intercept (and compute any intervals on that scale), or fit whatever more complicated model you like. To convert to a relationship between $$J$$ and $$\rho$$ you multiply each term in the $$v$$-equation by $$\rho$$ (including the intervals around the line).

If you fit a straight line to $$v$$ vs $$\rho$$ you get a quadratic through the origin for $$J$$ vs $$\rho$$... but the observations are weighted differently so the parameter estimates will be slightly different (an unweighted fit to $$J$$ vs $$\rho$$ overweights the points with large $$\rho$$ relative to the rest when they're relatively noisy/uncertain).

A caveat: all of the above discussion assumes that $$\rho$$ is essentially measured without noise. While it seems like it's going to be substantially less noisy/uncertain than measuring $$v$$ or $$J$$, I suspect it probably has some error and strictly speaking we should probably be looking at errors-in-variables or Model II regression (or at least to using fitting methods that are reasonably robust to x-errors). From the look of it I don't think anyone in the area is worrying about it, but I would at least be thinking about that issue and its potential impact in terms of bias and any effect on inference.

Edit 7:

For the quadratic fit to $$J$$, we can compute the empirical interval that contains 90% of the data values:

> quantile(abs(lm(J~rho+I(rho^2)-1)$residuals),0.9) 90% 0.2088087 Or we can compute one based on a normal approximation: > summary(lm(J~rho+I(rho^2)-1))$s*1.645
 0.2050971 (Or we could as easily do a $$t$$ interval)

But if we plot residuals it's obvious this interval isn't much use (which is why the normal approximation above looks pretty poor). We can see 90% (44 of 49) residuals are inside the bounds ... but for $$\rho$$ below 1, the residuals are way inside. The deviation around the model is obviously not constant as $$\rho$$ increases: On the other hand, if we fit a line to $$v$$ vs $$\rho$$, the intervals are more reasonable across the range of $$\rho$$, and we can calculate intervals in similar fashion to the ones above and multiply them by $$\rho$$ to get intervals that cover the same fraction of $$J$$-values overall but more accurately describe how far they tend to be from the curve.

> quantile(abs(residuals(lm(v~rho))),.9)
90%
0.1509835


Hello! Looks like that might be where they got the 0.15 from.

This is how that looks: Another way to generate an appropriate interval would be to model $$\log{(J)}$$ as a function of $$\rho$$; again, the residuals would have nearly constant spread.

--

Normally, of course, we're not so much interested in the fraction of the residuals that the model covers, but the coverage of an additional point (a prediction interval) - and that will be wider, because it has to take into account the uncertainty in the estimate of the curve.

• Glen, thank you. So you say we just fit a polynomial, exponential or a power trend line? (and that's what the article I showed has done?) – wmac Mar 29 '13 at 4:02
• I really don't know what they have done for sure. See my most recent edits. – Glen_b -Reinstate Monica Mar 29 '13 at 5:50
• I am sorry Glen. I think It has taken a lot of your time. It is very much appreciated. – wmac Mar 29 '13 at 12:18
• I think the best bet is to talk to people in the area of pedestrian flows and traffic modelling - especially the paper's authors. The models for this sort of data is well outside my expertise. – Glen_b -Reinstate Monica Mar 29 '13 at 12:22
• Basically I do a multi-agent simulation (i.e. agents which represent the pedestrian in crowd.) I needed to fit a line and then compare my line with those empirical data on ped-net website. I guess both red and green lines are good fits. But the red line is more realistic from my experience (since in high density crowd speed falls down and the flow will also gradually come down). At a point in time (around 5-6 peds/m^2) the crowd stalls. – wmac Mar 29 '13 at 12:27