I appears that bySwarm is right. Here are the results:
along the X axis is the total bonus over the six ability scores. Along the Y axis, the probability, obtained from 1 million runs. Results below a total bonus of +3 have been purged from the count, so the grand total of runs is less than the original 1 million.
It appears that the twelve 3d6 statistically produces a better total bonus than the 4d6 method.
This is the code to run the 4d6 case (in Python)
import sys
import random
count = {}
for i in xrange(1,1000000):
collection = []
for j in xrange(0,6):
extraction = [random.randint(1,6), random.randint(1,6), random.randint(1,6), random.randint(1,6) ]
#print extraction
collection.append( sum( sorted( extraction )[1:] ) )
#print collection
bonuses = map(lambda x: (x-10)/2, collection)
#print bonuses
total_bonus = sum(bonuses)
#print total_bonus
if total_bonus < 3:
#print "too low, excluded"
continue
if not count.has_key(total_bonus):
count[total_bonus]=0
count[total_bonus] += 1
total_extractions = sum(count.values())
for bonus,occurrences in sorted(count.items()):
print bonus,occurrences/float(total_extractions)
This is for the twelve 3d6 case:
import random
count = {}
for i in xrange(1,1000000):
collection = []
for j in xrange(0,12):
extraction = [random.randint(1,6), random.randint(1,6), random.randint(1,6) ]
collection.append( sum( extraction ) )
##print collection
collection = sorted(collection)[6:]
#print collection
bonuses = map(lambda x: (x-10)/2, collection)
#print bonuses
total_bonus = sum(bonuses)
#print total_bonus
if total_bonus < 3:
#print "too low, excluded"
continue
if not count.has_key(total_bonus):
count[total_bonus]=0
count[total_bonus] += 1
total_extractions = sum(count.values())
for bonus,occurrences in sorted(count.items()):
print bonus,occurrences/float(total_extractions)