Just some thoughts:
The script you linked allows a random X or Y offset up to 10m but doesn't check if the absolute offset is greater than 10m. For example, a Y offset of 9m with a X offset of 8m gives an offset of ~12m. I used to have VBA code that then looped the resulting offset to check if it was greater than the max dist and re-ran the randomizer but VBA is dead so that would be a waste of your time. You'd have to code the same kind of error-checker into the python script.
I do something similar but with client data so for you I'd suggest:
-Take your geocoded patient point file and add the centroid coordinates from the census block (join the census block data to the patient layer or spatial join if you geocoded to addresses)
-Make a XY event from the table of the geocoded based on the centroid values (to keep your source data intact)
-Add the randomizer script to model builder and run the XY even through the script
-Spatial join the randomized points to the census blocks
-Make two outputs- Matches to Census Blocks, Do not Match
-export the records that match the census block to an intermediate dataset
-take the remaining incorrectly placed patient points and re-run them through the randomizer step
-Spatial join again to export out the correctly placed points
-Then repeat the loop until the incorrect count is zero
You should be able to do all that in Model builder without coding (but don't quote me on that).
I can't help with the Python... I'm loyal to VBA
Have you been able to find any simpler approach? I am using claims data and have the same problem. I wonder whether overlapping points might be a problem for GWR.