|Nov 30th, 2012, 1:00:01 AM||#26|
Join Date: Feb 2012
Also GS, is this thread being continued?
Check out my RMT (Baton) Passport to Victory (Peaked #2)!
Want to read a warstory? Try this one: Dragons Dinos, and Donkeys (I mean ponies): an OU Warstory
|Jan 1st, 2013, 10:04:24 PM||#27|
Join Date: Jun 2012
Hi all - I'm an Australian electrical engineering student who has been mostly lurking on Smogon for the past 6 months or so. I'm quite passionate about interesting mathematical problems though, so this thread caught my attention.
Anyway, I believe I have an analytical solution to Problem 1 (in its entirety). I've used MATLAB extensively in developing the solution, so pretty much all of the calculations are done via MATLAB code. While the maths behind said code is generally pretty simple - assuming you're familiar with Markov chains, matrices in general and basic combinatorics, you're probably alright - it gets a bit messy along the way so bear with me. Because the solution requires dealing with large matrices, I've been unable to provide a general solution, only a specific one dealing with the 7 stat, 13 state Moody specific problem. As an engineering student, that doesn't bother me so much, but I know the more pure maths types out there might be bothered by it so I figured I'd include a disclaimer :P On the note of disclaimers, I'll include a tl;dr warning now - this is not a short post.
As quite correctly stated by alkinesthetase earlier in the thread, to deal with a pure 7 variable, 13 state problem you'd be looking at ~62-63 million states (and thus a ~4*10^15 element transition matrix). Even MATLAB won't happily deal with matrices in that order of magnitude, and even if it could the processing time would be horrendous (several hours minimum to just generate a list of a list of all possible states, let alone fill in that transition matrix). This is where one of only about two major mathematical leaps required to solve this problem practically steps in. Based on the definition of the problem, all stats are indistinguishable. As such, a stat spread where after the first turn Attack is raised to +2 and Defence is dropped to -1 is considered equivalent to Evasion rising to +2 and Speed to -1.
With this in mind, we can take a slightly different approach to the solution. Instead of considering all 13^7 stat modifier combinations, we can consider only those which are unique under the conditions mentioned above. A brute force way to achieve this would be to generate all 13^7 possible combinations and then filter those with duplicates, but as mentioned earlier just generating those states takes a prohibitively long amount of time. The code I used in the solution stores the states in a 50388x7 matrix, with each column containing one arbitrary stat value. Instead of working with -6,-5,...,+5,+6 I shifted the 13 states to 1-13 (I chose 1 as the minimum value because of the way array indexing works in MATLAB, with the first element listed as 1 instead of 0 as in some other programming languages). As such, each row of the matrix contains 7 values in the range 1 to 13.
The code functions by counting up from left to right (i.e. begin with all values at 1, increment the leftmost value until it reaches 13, then overflow into the second column stat). However, in order to eliminate non-unique modifier combinations, after an overflow each stat is set to a minimum value of the stat immediately to its right (essentially forcing descending order of the modification values). I can't offer a rigorous mathematical proof of why this will produce only the unique, order independent combinations. However, hopefully you can see intuitively why it is the case - for example, after 13 1 1 1 1 1 1 overflows the next combination listed is 2 2 1 1 1 1 1, since 1 2 1 1 1 1 1 is equivalent to 2 1 1 1 1 1 1, which would already have been generated previously. The code written for this is shown in the spoiler below.
I mentioned the magic number 50388 earlier without much context. The output of the code above is indeed a matrix of that size (MATLAB automatically increases the size of a matrix if it's too small, so unlike other languages selecting an incorrect initial size for the variable doesn't automatically result in an incorrect code function, just somewhat slower code). However, there are a few ways to verify that that is the required number of combinations. One of these is covered later, while developing a formula for the location of any given combination in the list of all combinations. There is an alternative that I stumbled across while first approaching this problem - I've spoilered it below because it isn't actually essential to the overall solution.
With the list of states generated, I could now attempt to generate the 50388x50388 transition matrix required to solve the problem as a Markov chain. I pursued a number of ultimately inefficient paths to solve this part of the problem, mostly because of the complexity I perceived to be involved with probably the single most interesting aspect of the overall Problem 1 solution. The intuitively simple way to approach this step is to look at each state, and generate all possible states to which it can transition by increasing one stat +2 and decreasing another by -1.
There are 3 exceptional cases which needed to be handled - firstly, the possibility that all stats were at 13, in which case there was a probability of 1 that it would transition to 13 13 13 13 13 13 12. Similarly, all stats at 1 had to transition to 3 1 1 1 1 1 1. The case with 6 stats at 1 and 1 stat at something other than 1 was a little more interesting. I made the assumption when solving this problem that, at a code level, Moody selects that stat that it is going to raise before it selects the stat it will decrease. By definition, one has to occur first, since true simultaneity can't exist if they are prevented from selecting the same stat. As such, while most of the time for this case one of the stats at 1 will increase to 3 and the non-1 stat decreased, there is a 1/7 chance that the increase will select the non-1 stat, with no stat then able to be decreased.
For all cases other than the exceptions, the code was relatively simple. Check all possible combinations of indices, if the two indices are not the same, the index chosen for the increasing stat doesn't point to a stat at 13 and the index chosen for the decreasing stat doesn't point to a stat at 1, increase and decrease the respective stats and store the new state (reordered to be descending again) in a temporary matrix. Once this is done, check the matrix for duplicates, clear out the duplicates and store an additional vector containing the number of duplicates each state had. Essentially, this gives me a matrix containing all the possible outcome states from the given initial state, as well as a vector containing the relative probability that a given outcome state will be selected. I've left the code for this until the end of this section so that it can be seen in its entirety.
This leads to the single most interesting part of the problem, and by far the most interesting in my opinion - how do I place those probabilities into the transition matrix? While I have the outcome states listed, I need to convert them into a matrix index before I can actually store the probabilities. I foresaw this problem early on, which is largely why I tried so many other methods before tackling this one. A brute force lookup table is incredibly slow given the number of possible states, so that was ruled out rather quickly. However, there is a way.
Going back to the process by which the state list matrix was filled, I was able to observe a pattern. For the first stat alone, each increase corresponded to 1 new state - simple enough. The second stat was slightly more complex. If the second stat was at 1, there were a complete 13 states the first stat could occupy. This decreased to 12 for the second stat at 2 because the first stat would then start counting at 2, dropping to 11 for 3, etc. Essentially, the total number of states up to 13 13 1 1 1 1 1 was then given by the 13th triangular number. Increasing up to 13 13 2 1 1 1 1 added an additional number of states equal to the 12th triangular number, with a similar decay occurring for each increase to that of the second stat. In the end, the total number of states for the third stat was the sum of the sum of the integers 1-13. Following the pattern through, the fourth stat was the sum of the sum of the sum, the fifth stat the fourth tier sum, through to the seventh stat at the sixth tier sum. Sure enough, running that calculation - the sixth tier sum from 1-13 - once again generates that number 50388.
To me, this was an incredibly cool mathematical pattern, and even if I hadn't managed to convert it into an overall solution would have made the entire process worthwhile. It required a little additional tweaking to be entirely useful, with several separate calculations done to include the account for intermediate states between 13...13 (x) 1...1 states. In isolation, this code is shown below.
With this code written, it became rather elementary to store the probabilities. The probability of a given transition was equal to the relevant entry in the duplicates vector divided by the sum of the vector (the relative probability over the sum of all relative probabilities). The code above could then be run on the relevant member of the outcome state matrix to give the transition matrix index in which the probability would be stored.
The overall transition matrix code is shown below, including a duplication of the code mentioned previously in the correct location. In developing this code, I've assumed a column probability vector and pre-multiplying transition matrix since that's how I was taught way back in my high school days - this means that the columns, not the rows, of the transition matrix will sum to 1. Transition matrix columns correspond to the initial state; rows to the outcome state. It's rather important to note that, even with the dramatic size reduction achieved by this method of solution, the transition matrix would still be too large to store normally in MATLAB. As such, it is implemented in the code specifically as a sparse matrix, which fits well given how few non-zero elements it has. This does become important in subsequent parts, however. MATLAB's coding language allows sparse matrices to be treated essentially the same as normal matrices, however, so it doesn't actually effect how this particular section of the code is written. I should add that I'm no coding expert - I'd never coded in MATLAB prior to this, so I've no doubt that there are inefficiencies that could be further improved upon. However, shown below is the code that I used - it took around 10 minutes to run, which is reasonable enough given the magnitude of the matrices involved in my opinion.
With the list of states and the transition matrix developed, it became elementary to solve the remainder of the problem. The one important choice was between implementing the transitions as a power of a matrix or as an iterative procedure. Raising the transition matrix to a large power results in most of the zero elements becoming non-zero, removing the benefits of the sparse matrix definition that allowed MATLAB to handle the matrix in the first place. On the other hand, by iterating the process, the matrix always remained sparse, with only the probability vector becoming predominantly non-zero.
At each iteration, for the part a) solution the probability vector needed to be converted to a set of probabilities for an arbitrary stat being at a particular modification level, not at a particular state. This was relatively simple - firstly, an output matrix 51x13 was defined (51 to allow for 50 turns plus the initial state; 13 to allow for each possible modification level). At each iteration, for each state, the probability of being at that particular state was multiplied by the probability of an arbitrary stat from that state being at a given modification level (expressed in a matrix), with the result added to the relevant output cell. The code is once again shown in the spoiler below.
As for results, the full probability distribution is shown in the spoiler below.
However, I will highlight this particular section, the results specifically for turn 50:
-6 0.000002592372339 -5 0.000010761952707 -4 0.000046592692124 -3 0.000146708498645 -2 0.000420155868716 -1 0.001042149121095 +0 0.002386656940549 +1 0.005145830941225 +2 0.011225952420392 +3 0.026896542574243 +4 0.076036579426630 +5 0.252086671245168 +6 0.624552805946163
Part b) was handled with very similar code to part a) - I calculated them as part of the same MATLAB script. Once again, iterations were used to calculate the probability vector for each stage. However, this time probability of a given state needed to be converted to the probability of a particular number of positive boosts. Once again, these were expressed in a matrix and multiplied by the state probability vector, giving a 51x43 matrix (50 turns plus initial state; 7 stats at a maximum of +6 plus the 0 state). An additional column was then added for the average number of positive boosts, plus a column for the average Base Power of the move Stored Power on that turn. The code and total probability distribution are shown below.
Once again, highlighting only the turn 50 results, the data agrees well with ganj4lF's simulations.
Turn +'ve Boosts SP BP 0 0 20.0000 1.0000 2.0000 60.0000 2.0000 3.7143 94.2857 3.0000 5.2041 124.0816 4.0000 6.5287 150.5736 5.0000 7.7348 174.6951 6.0000 8.8561 197.1211 7.0000 9.9155 218.3105 8.0000 10.9280 238.5593 9.0000 11.9027 258.0531 10.0000 12.8454 276.9071 11.0000 13.7598 295.1958 12.0000 14.6486 312.9711 13.0000 15.5136 330.2724 14.0000 16.3566 347.1326 15.0000 17.1790 363.5806 16.0000 17.9821 379.6425 17.0000 18.7671 395.3419 18.0000 19.5350 410.7005 19.0000 20.2869 425.7381 20.0000 21.0236 440.4727 21.0000 21.7460 454.9208 22.0000 22.4549 469.0976 23.0000 23.1508 483.0166 24.0000 23.8345 496.6903 25.0000 24.5065 510.1295 26.0000 25.1672 523.3443 27.0000 25.8172 536.3435 28.0000 26.4568 549.1352 29.0000 27.0863 561.7265 30.0000 27.7062 574.1234 31.0000 28.3165 586.3309 32.0000 28.9176 598.3518 33.0000 29.5094 610.1874 34.0000 30.0918 621.8369 35.0000 30.6649 633.2983 36.0000 31.2284 644.5688 37.0000 31.7823 655.6462 38.0000 32.3265 666.5297 39.0000 32.8610 677.2200 40.0000 33.3859 687.7181 41.0000 33.9012 698.0234 42.0000 34.4066 708.1312 43.0000 34.9015 718.0295 44.0000 35.3849 727.6981 45.0000 35.8554 737.1076 46.0000 36.3110 746.2205 47.0000 36.7497 754.9937 48.0000 37.1691 763.3814 49.0000 37.5669 771.3389 50.0000 37.9413 778.8257
Firstly, to anyone who made it this far through the wall of text above, thanks for reading. Secondly, thanks to Great Sage for posing such an interesting problem, alkinesthetase for having a shot at it and providing his time to let me bounce ideas around and ganj4lF for providing the aforementioned simulation data. As I said at the beginning, I love these sorts of problems, so even if the outcome is entirely pointless it was fun just exploring the maths involved with something like Moody.
I also have a set of data for a hypothetical (and utterly ridiculous) situation in which someone uses Baton Pass to bring a Moody Pokemon in with -6 in all stats. I'll probably post it at some point in the future; however, I've been talking to Great Sage about some inconsistencies with it (my guess, most probably associated with my assumption about the increase stat always being selected first), so I'll wait until that's all cleared up before posting.
Anyway, I hope this has been of interest to someone - back to lurking for me :D
Last edited by Tsaeb XIII; Jan 1st, 2013 at 11:47:08 PM.