Here's a simple program which tells us what the distribution of each individual stat (from highest to lowest) in your array is:
ARRAY: 7d[highest 3 of 4d{3..6}]
output 1@ARRAY
output 2@ARRAY
output 3@ARRAY
output 4@ARRAY
output 5@ARRAY
output 6@ARRAY
In order to generate a sequence of 7 rolled ability scores, we can simply use the function which generates one ability score and roll it as a die 7 times - hence 7d[highest 3 of 4d{3..6}]
. Then, we can inspect the first six elements of that sequence to get the six ability scores you will keep, since sequences are by default sorted in descending order. The 7th score (7@ARRAY
) is the lowest score and so we ignore (discard) it.
Here's a program which compares results from simply rolling 6 times 3d6 to your method of determining stats. As given it only shows the comparison of the highest score from each distribution, but I trust altering it to compare the other results is obvious enough (I have left them out to avoid clutter):
ORIGINAL: 6d(3d6)
NEW: 7d[highest 3 of 4d{3..6}]
loop N over {1} {
output N@ORIGINAL named "[N]@3d6"
output N@NEW named "[N]@7d(4d6r1,2)"
}

That gives us a graph outlining the expected distribution for each ability score, but the distributions are separate. If you want to aggregate those into a single distribution giving you the expected distribution of any one randomly selected score from your kept array, we need to get a bit funky with a function, as in this program:
function: select INDEX:n DISTRIBUTION:d {
result: INDEX@DISTRIBUTION
}
ARRAY: 7d[highest 3 of 4d{3..6}]
output [select 1d6 ARRAY]
In this case we use a function in order to randomly select an index from the array. We have to do it this way because anydice expects an index to a sequence to be a number or another sequence (where it would sum the selected elements), so it balks if we try to directly index the sequence with a die by doing something like (1d6)@ARRAY
. However, if we use a function which casts a die roll input to a number (INDEX:n
), we can work around this limitation. The output in this case gives us the probability distribution of results we'd see if we randomly selected one of the first six ability scores in the sequence we generated.
Using this program, here's a graph which compares that to other methods of generating ability scores:

So we can see, for instance, that adding the extra step of rolling an array of seven scores and dropping the lowest bumps your average result up from 14.62 to 15.06.