A little but time consuming experiment first.  You'll not need to do it, as it is only for finding a model.

Let's see how the Min value of the plot intensity varies with Opacity:

    theData = RandomReal[NormalDistribution[], {2000, 2}];
    f[x_] := f[x] = 
       Min[Norm /@ Flatten[ImageData@Rasterize[
            ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200, 
                       PlotStyle -> {Black, Opacity[x]}, Axes -> False]], 1]];
    Plot[f[x] , {x, 0, .4}, PlotRange -> Full]

![Mathematica graphics](http://i.stack.imgur.com/rP4Tc.png)

So, it is an exponential. (**Note**: in the last edit to this post I got rid of the exponential model, fitting a few points with an Interpolation, and it works nice)

Let's fit it:

    data = Table[{i, f[i]}, {i, 0, 1, .1}]

    model = a Exp[b x];
    fit = FindFit[data, model, {a, b}, x];
    modelf = Function[{t}, Evaluate[model /. fit]]

    Show[ListPlot@data, Plot[modelf[x], {x, 0, 1}]]

![Mathematica graphics](http://i.stack.imgur.com/vEOEz.png)


Now you are ready to set the min value of the brightness of the plot to whatever you want:

(The Sqrt@3 is a normalization factor for the intensity of the {1,1,1} RGB pixel.)

Let's use it:

    opac = x /. Solve[# == a E^(b x)/Sqrt@3, x] /. fit & /@ {1/2, 1/4, 1/20, 1/200}
    
    ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200, 
       PlotStyle -> {Black, Opacity[#[[1]]]}, Axes -> False] & /@ opac

![Mathematica graphics](http://i.stack.imgur.com/3cIsR.png)

**Edit**

Let's pack that all together in a function and plot two very different point sets with the same maximum darkness.

    ClearAll[opa];
    Options[opa] = Options[ListPlot];
    opa[desiredOpacity_, points_, opts : OptionsPattern[]] :=
     Module[{f, a, b, model, fit, modelf, x},
      f[x_] := 
       f[x] = Min[
         Norm /@ Flatten[
           ImageData@
            Rasterize[
             ListPlot[points, Axes -> False, 
              PlotStyle -> {Black, Opacity[x]}]], 1]];
      
      model = a Exp[b x];
      fit = FindFit[Table[{i, f[i]}, {i, 0, 1, .1}], model, {a, b}, x];
      modelf = Function[{t}, Evaluate[model /. fit]];
      Return[x /. Quiet@Solve[# == modelf[x]/Sqrt@3, x][[1]] /. fit &@
        desiredOpacity];
      ]
    
    theData  = RandomReal[NormalDistribution[], {2000, 2}];
    theData1 = RandomReal[NormalDistribution[], {10000, 2}];
    
    opad = opa[.5,  theData,  AspectRatio -> Automatic, ImageSize -> 200];
    opad1 = opa[.5, theData1, AspectRatio -> Automatic, ImageSize -> 200];
    
    Grid[{{ListPlot[theData, Axes -> False, PlotStyle -> {Black, Opacity[opad]},
                             AspectRatio -> Automatic, ImageSize -> 200], 
           ListPlot[theData1,Axes -> False, PlotStyle -> {Black, Opacity[opad1]}, 
                             AspectRatio -> Automatic, ImageSize -> 200]}}, 
       Frame -> All]

![Mathematica graphics](http://i.stack.imgur.com/6TtVg.png)


The same, but darker 

![Mathematica graphics](http://i.stack.imgur.com/qEtLT.png)


Remember that the default plot is:


![Mathematica graphics](http://i.stack.imgur.com/Atv3d.png)


**Edit**

Answering @Oleksandr comments, the following does not assume an exponential model:

    ClearAll[opa];
    Options[opa] = Options[ListPlot];
    opa[desiredOpacity_, points_, opts : OptionsPattern[]] := 
     Module[{f, model, x}, 
      f[x_] := f[x] = 
        Min[Norm /@ Flatten[ImageData@ Rasterize[
             ListPlot[points, Axes -> False, PlotStyle -> {Black, Opacity[x]}]], 1]];
      model = Interpolation[Table[{i, f[i]}, {i, 0, 1, .1}]];
      x /. FindRoot[desiredOpacity == model[x]/Sqrt@3, {x, 0, 1}][[1]]]