Nowhere in my previous post am I assuming PPLNS. I am actually assuming straight PPS. Please read it and try to understand it. Or at least look at the python script.
You were assuming straight PPS, I was assuming per block proportional. Do you see? We will never convince each other because we are both right, in the respective systems we think we are operating under.
Sorry, but you are wrong. You were not simulating per block proportional, you were assuming a wrong distribution for the time between shares.
Thank you for the detailed, thorough response. Give me a chance to digest everything you said. You are clearly more advanced in mathematics and statistics than I (i mean that genuinely not sarcastically) and I would love to learn from you. However, in the mean time, let me ask you this.
If you can't precalculate (you say precalculate, I say predict) the time is takes to solve a block, how is bitcoin (or whatever coin) doing it? Bitcoin itself has to predict the time it will take the network to solve a block, to know how to adjust the difficulty. It MUST maintain an average block find time of 10 minutes, or quite literally the entire thing will come crashing to a screeching halt.
If bitcoin couldn't accurately predict the time it takes to solve a block given X hashpower, we wouldn't be having this conversation because bitcoin would have died, and none of the alt coins that we love to mine at middlecoin would have ever existed.
I elaborated a bit at the end, because that's what I thought was most likely to be read:
Now on to your simulation: You are using the wrong distribution. If you don't know why, read my post again from the top, read the wiki page or watch the video I linked above.
The probability distribution of the times between finding shares (again, this is a poisson process) is simply e^(-t/T), where T is the average time between shares. Which has exactly the property that you can just restart it at all points without changing anything.
edit to clarify: the chance P of the time you need to solve a block is P(t) = e^(-t/T). For generating t, you can for example use numpy.random.exponential(scale=T)