I think the main problem is one of dynamic range, your algorithm is probably right but you're working on the wrong type of data.

A point light source that would otherwise clip and go pure white gets spread over a larger area by a defocussed lens, so that it forms a disc that isn't as bright and therefore doesn't clip.

That's why you get those nice circles in your real bokeh image. If you clip the signal (making it less bright than it otherwise would be and then spread it out with your bokeh simulation you get a dim circle (or hexagon, or whatever) that doesn't stand out and thus doesn't look realistic.

What you have in a real image chain is:

bokeh (from the lens) -> digitisation (clipping) -> gamma correction & dynamic range compression

What you are doing is

sharp image -> digitisation (clipping) -> gamma correction & dynamic range compression -> bokeh simulation

You wont get the correct result because you're not working with linear data.

What you can do is attempt to linearise the data, replace any dynamic range that has been lost to clipping, perform your bokeh simulation, and then redo the nonlinear operations!

Here's an example. I've started with an HDR image that has been tonemapped, giving a highly nonlinear result. This is the worst type of image to attempt bokeh simulation with!

Doing a standard convolution operation to simulate bokeh (using photoshop's lens blur tool) yields this result, which is very similar to what you are getting:

To get a better result, I applied an extreme curve to try and get the image back to roughly what it would have been before tonemapping, where the highlights are much, much brighter than the rest of the image. I did this by with the levels tool, pushing the centre input a long way to the right, from 1.0 to about 0.2). I then applied the lens blur tool, just as before. Finally I applied an extreme curve in the opposite direction to the first curve. The result, whilst a long way from perfect, looks much more like real lens bokeh:

 

If you're doing this in code, then try cubing each value, then applying your bokeh simulation routine, then take the cube root of each value. You should see an improvement. It might take some tweaking.

tl;dr even if you have implemented a perfect mathematical model of bokeh, it must be applied on unclipped linear data. If you apply the same calculations to heavily modified data (even a standard in camera JPEG is heavily modified from a mathematical point of view) you will get a very different result.