Some mathematical problems

[Kidogo]

Originally Posted by Fat Antar
I solved Problem 3 months ago.
Only 3 months ago? We solved this 3 years ago.

Hey, my answer was literally one away from the answer here--I probably made an OBOE error. It supports my method though!

Also GS, is this thread being continued?

[Tsaeb XIII]

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.


Foundation Calculations:

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.

...

---

% Define stat modification variables
s1=1;
s2=1;
s3=1;
s4=1;
s5=1;
s6=1;
s7=1;
% Define counting variables
k=1; % Output matrix index variable
% Define output matrix
State_list=zeros(50388,7);
while s7<14
s6=s7;
while s6<14
s5=s6;
while s5<14
    s4=s5;
    while s4<14
        s3=s4;
        while s3<14
            s2=s3;
            while s2<14
                s1=s2;
                while s1<14                      
                    State_list(k,:)=[s1 s2 s3 s4 s5 s6 s7];                              
                    k=k+1;                                                       
                    s1=s1+1;
                end
                s2=s2+1;
            end
            s3=s3+1;
        end
        s4=s4+1;
    end
    s5=s5+1;
end
s6=s6+1;
end
s7=s7+1;
end

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.

...

---

When I first decided to compress non-unique states into a single unique state, I was using a different way of storing my results, with 13 columns containing the number of stats at each boost level (i.e. 13 values from 0 to 7 that sum to 7). This proved to be an entirely unwieldy way of actually solving the problem, but it did offer a useful way of calculating how many distinct combinations there were. Because there is a finite number of ways to break 7 up into a sum of values from 1 to 7 (Google 'partitions of 7' and look for the mathematics related results), you can essentially break the problem into a bunch of small combinatorics calculations. In each case, nCr/Pascal's triangle gives the number of combinations if there weren't duplicates, and dividing by the factorial of the number of duplicates in each case gives the true number of combinations. Adding the combinations from each case once again gives the magic number 50388.

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.

...

---

 Sb=0;
 S1=0;
 S2=0;
 S3=0;
 S4=0;
 S5=0;
 S6=0;
 S7=0;
 State_index=0;
 
 initial_state=outcome_states(a,:); % using initial_state as a temp variable
 
 case_no=0;
 max_flag=0;
 n=7;
 while n>0
     if initial_state(1,n)==1
         case_no=n;
     end
     n=n-1;
 end
 if sum(initial_state==13)==7
     case_no=0;
 elseif case_no==0
     case_no=8;
 end
 case_no=case_no-1;
 
 % Note that the commented summing loops are those that would
 % produce the Sb values given below, but are not required for the
 % function
 switch case_no
     case 0
         Sb=1;
 
     case 1
         Sb=1;
 
     case 2
         Sb=13;
 
     case 3
%               k2=13;
%               while k2>0
%                   k1=k2;
%                   while k1>0
%                       Sb=Sb+1;
%                       k1=k1-1;
%                   end
%                   k2=k2-1;
%               end
             Sb=91;
 
     case 4
%               k3=13;
%               while k3>0
%                   k2=k3;
%                   while k2>0
%                       k1=k2;
%                       while k1>0
%                           Sb=Sb+1;
%                           k1=k1-1;
%                       end
%                       k2=k2-1;
%                   end
%                   k3=k3-1;
%               end
             Sb=455;
 
     case 5
%               k4=13;
%               while k4>0
%                   k3=k4;
%                   while k3>0
%                       k2=k3;
%                       while k2>0
%                           k1=k2;
%                           while k1>0
%                               Sb=Sb+1;
%                               k1=k1-1;
%                           end
%                           k2=k2-1;
%                       end
%                       k3=k3-1;
%                   end
%                   k4=k4-1;
%               end
             Sb=1820;
 
     case 6
%               k5=13;
%               while k5>0
%                   k4=k5;
%                   while k4>0
%                       k3=k4;
%                       while k3>0
%                           k2=k3;
%                           while k2>0
%                               k1=k2;
%                               while k1>0
%                                   Sb=Sb+1;
%                                   k1=k1-1;
%                               end
%                               k2=k2-1;
%                           end
%                           k3=k3-1;
%                       end
%                       k4=k4-1;
%                   end
%                   k5=k5-1;
%               end
             Sb=6188;
 
     case 7
%               k6=13;
%               while k6>0
%                   k5=k6;
%                   while k5>0
%                       k4=k5;
%                       while k4>0
%                           k3=k4;
%                           while k3>0
%                               k2=k3;
%                               while k2>0
%                                   k1=k2;
%                                   while k1>0
%                                       Sb=Sb+1;
%                                       k1=k1-1;
%                                   end
%                                   k2=k2-1;
%                               end
%                               k3=k3-1;
%                           end
%                           k4=k4-1;
%                       end
%                       k5=k5-1;
%                   end
%                   k6=k6-1;
%               end
             Sb=18564;
 
     case -1
     max_flag=1;
     k7=13;
%               while k7>0;
%                   k6=k7;
%                   while k6>0;
%                       k5=k6;
%                       while k5>0;
%                           k4=k5;
%                           while k4>0;
%                               k3=k4;
%                               while k3>0;
%                                   k2=k3;
%                                   while k2>0;
%                                       k1=k2;
%                                       while k1>0;
%                                           Sb=Sb+1;
%                                           k1=k1-1;
%                                       end
%                                       k2=k2-1;
%                                   end
%                                   k3=k3-1;
%                               end
%                               k4=k4-1;
%                           end
%                           k5=k5-1;
%                       end
%                       k6=k6-1;
%                   end
%                   k7=k7-1;
%               end
             Sb=50388; 
 end
 
 if case_no>=7
     k7=12;
     while k7>(14-initial_state(1,7))
         k6=k7;
         while k6>0
             k5=k6;
             while k5>0
                 k4=k5;
                 while k4>0
                     k3=k4;
                     while k3>0
                         k2=k3;
                         while k2>0
                             k1=k2;
                             while k1>0
                                 S7=S7+1;
                                 k1=k1-1;
                             end
                             k2=k2-1;
                         end
                         k3=k3-1;
                     end
                     k4=k4-1;
                 end
                 k5=k5-1;
             end
             k6=k6-1;
         end
         k7=k7-1;
     end
 end
 if case_no>=6
     k6=14-initial_state(1,7)-(initial_state(1,7)==1);
     while k6>(14-initial_state(1,6))
         k5=k6;
         while k5>0
             k4=k5;
             while k4>0
                 k3=k4;
                 while k3>0
                     k2=k3;
                     while k2>0
                         k1=k2;
                         while k1>0
                             S6=S6+1;
                             k1=k1-1;
                         end
                         k2=k2-1;
                     end
                     k3=k3-1;
                 end
                 k4=k4-1;
             end
             k5=k5-1;
         end
         k6=k6-1;
     end
 end
 if case_no>=5
     k5=14-initial_state(1,6)-(initial_state(1,6)==1);
     while k5>(14-initial_state(1,5))
         k4=k5;
         while k4>0
             k3=k4;
             while k3>0
                 k2=k3;
                 while k2>0
                     k1=k2;
                     while k1>0
                         S5=S5+1;
                         k1=k1-1;
                     end
                     k2=k2-1;
                 end
                 k3=k3-1;
             end
             k4=k4-1;
         end
         k5=k5-1;
     end
 end
 if case_no>=4
     k4=14-initial_state(1,5)-(initial_state(1,5)==1);
     while k4>(14-initial_state(1,4))
         k3=k4;
         while k3>0
             k2=k3;
             while k2>0
                 k1=k2;
                 while k1>0
                     S4=S4+1;
                     k1=k1-1;
                 end
                 k2=k2-1;
             end
             k3=k3-1;
         end
         k4=k4-1;
     end
 end
 if case_no>=3
     k3=14-initial_state(1,4)-(initial_state(1,4)==1);
     while k3>(14-initial_state(1,3))
         k2=k3;
         while k2>0
             k1=k2;
             while k1>0
                 S3=S3+1;
                 k1=k1-1;
             end
             k2=k2-1;
         end
         k3=k3-1;
     end
 end
 if case_no>=2
     k2=14-initial_state(1,3)-(initial_state(1,3)==1);
     while k2>(14-initial_state(1,2))
         k1=k2;
         while k1>0
             S2=S2+1;
             k1=k1-1;
         end
         k2=k2-1;
     end
 end
 if case_no>=1
     S1=initial_state(1,1)-initial_state(1,2)+(max_flag==0)-(case_no==1);
 end
 
 State_index=Sb+S7+S6+S5+S4+S3+S2+S1;

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.

...

---

% This script requires Permutation_Generator to be run first, generating
% the matrix State_list containing the possible states. Dimensions of State_list are
% 50388 x 13, with each row containing a distinct state.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variable Definitions:
% Define transition matrix
Transition_matrix=sparse(1,1,0,50388,50388);
% Define array variables
initial_state=zeros(1,7); % storage vector for current state
outcome_states=zeros(1,7); % storage matrix for possible outcome states
relative_probabilities=zeros(1); % vector for storing number of occurances of each state
% Define scalar variables
i=1; % Index variable for current state (i.e. column of transition matrix)
t=0; % Total number of possible state changes from the current state
inc=0; % Counting variables for 
dec=0; % generation  
j=0; % of
a=0; % output
b=0; % states
Sb=0; % Base (13's followed by 1's sum)
S1=0; % Increase on base due to s1
S2=0; % Increase on base due to s2
S3=0; % Increase on base due to s3
S4=0; % Increase on base due to s4
S5=0; % Increase on base due to s5
S6=0; % Increase on base due to s6
S7=0; % Increase on base due to s7
State_index=0; % Location of given state in list of all states
k1=0; % Counting
k2=0; % variables
k3=0; % for
k4=0; % state
k5=0; % location
k6=0; % in
k7=0; % list
n=0; % Counting variable for identifying case while locating state
case_no=0; % Case identification variable while locating state
max_flag=0; % Flag variable for all 13 case while locating state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main Code
while i<50389
 % Output state generation setup
 
 initial_state=State_list(i,:); % Retrieve current state
 outcome_states=zeros(1,7); % Clear outcome state matrix
 relative_probabilities=zeros(1);
 inc=1; % Clear
 dec=1; % counting
 j=1; % variables
 a=1;
 b=1;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Generate possible output states
 
 if sum(initial_state==13)==7
     outcome_states=[13 13 13 13 13 13 12];
     relative_probabilities=[1];
     j=2;
  elseif sum(initial_state==1)==6
     outcome_states(1,:)=fliplr(sort([1 1 1 1 1 3 (initial_state(1,1)-1)]));
     if initial_state(1,1)==12
         outcome_states(2,:)=[(initial_state(1,1)+1) 1 1 1 1 1 1];
         relative_probabilities=[(6/7);(1/7)];
         j=3;
     elseif initial_state(1,1)==13
         relative_probabilities=[1];
         j=2;
     else
         outcome_states(2,:)=[(initial_state(1,1)+2) 1 1 1 1 1 1];
         relative_probabilities=[(6/7);(1/7)];
         j=3;
     end
 elseif sum(initial_state==1)==7
     outcome_states=[3 1 1 1 1 1 1];
     relative_probabilities=[1];
     j=2;
 else
     while inc<8 % Stat to be increased
         dec=1;
         while dec<8 % Stat to be decreased
             if (inc~=dec)*(initial_state(1,inc)~=13)*(initial_state(1,dec)~=1)                
                 outcome_states(j,:)=initial_state;
 
                 if outcome_states(j,inc)==12
                     outcome_states(j,inc)=13;
                 else
                     outcome_states(j,inc)=outcome_states(j,inc)+2;
                 end
 
                 outcome_states(j,dec)=outcome_states(j,dec)-1;
 
                 outcome_states(j,:)=fliplr(sort(outcome_states(j,:)));
                 j=size(outcome_states,1)+1;
 
             end            
             dec=dec+1;
         end
         inc=inc+1;
     end
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Convert output states to unique options only
 
 a=1;
 while a<j
    if a==1
        relative_probabilities(1,1)=1;
    end
    b=a-1;
    while b>0
        if isequal(outcome_states(a,:),outcome_states(b,:))
            outcome_states(a,:)=[];
            j=j-1;
            relative_probabilities(b,1)=relative_probabilities(b,1)+1;
            a=a-1;
            break
        elseif b==1
            relative_probabilities(a,1)=1;
        end
        b=b-1;               
    end
    a=a+1;
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Place probabilities in transition matrix
 
 a=1;
 while a<j
     Sb=0;
     S1=0;
     S2=0;
     S3=0;
     S4=0;
     S5=0;
     S6=0;
     S7=0;
     State_index=0;
 
     initial_state=outcome_states(a,:); % using initial_state as a temp variable
 
     case_no=0;
     max_flag=0;
     n=7;
     while n>0
         if initial_state(1,n)==1
             case_no=n;
         end
         n=n-1;
     end
     if sum(initial_state==13)==7
         case_no=0;
     elseif case_no==0
         case_no=8;
     end
     case_no=case_no-1;
 
     % Note that the commented summing loops are those that would
     % produce the Sb values given below, but are not required for the
     % function
     switch case_no
         case 0
             Sb=1;
 
         case 1
             Sb=1;
 
         case 2
             Sb=13;
 
         case 3
%               k2=13;
%               while k2>0
%                   k1=k2;
%                   while k1>0
%                       Sb=Sb+1;
%                       k1=k1-1;
%                   end
%                   k2=k2-1;
%               end
                 Sb=91;
 
         case 4
%               k3=13;
%               while k3>0
%                   k2=k3;
%                   while k2>0
%                       k1=k2;
%                       while k1>0
%                           Sb=Sb+1;
%                           k1=k1-1;
%                       end
%                       k2=k2-1;
%                   end
%                   k3=k3-1;
%               end
                 Sb=455;
 
         case 5
%               k4=13;
%               while k4>0
%                   k3=k4;
%                   while k3>0
%                       k2=k3;
%                       while k2>0
%                           k1=k2;
%                           while k1>0
%                               Sb=Sb+1;
%                               k1=k1-1;
%                           end
%                           k2=k2-1;
%                       end
%                       k3=k3-1;
%                   end
%                   k4=k4-1;
%               end
                 Sb=1820;
 
         case 6
%               k5=13;
%               while k5>0
%                   k4=k5;
%                   while k4>0
%                       k3=k4;
%                       while k3>0
%                           k2=k3;
%                           while k2>0
%                               k1=k2;
%                               while k1>0
%                                   Sb=Sb+1;
%                                   k1=k1-1;
%                               end
%                               k2=k2-1;
%                           end
%                           k3=k3-1;
%                       end
%                       k4=k4-1;
%                   end
%                   k5=k5-1;
%               end
                 Sb=6188;
 
         case 7
%               k6=13;
%               while k6>0
%                   k5=k6;
%                   while k5>0
%                       k4=k5;
%                       while k4>0
%                           k3=k4;
%                           while k3>0
%                               k2=k3;
%                               while k2>0
%                                   k1=k2;
%                                   while k1>0
%                                       Sb=Sb+1;
%                                       k1=k1-1;
%                                   end
%                                   k2=k2-1;
%                               end
%                               k3=k3-1;
%                           end
%                           k4=k4-1;
%                       end
%                       k5=k5-1;
%                   end
%                   k6=k6-1;
%               end
                 Sb=18564;
 
         case -1
         max_flag=1;
         k7=13;
%               while k7>0;
%                   k6=k7;
%                   while k6>0;
%                       k5=k6;
%                       while k5>0;
%                           k4=k5;
%                           while k4>0;
%                               k3=k4;
%                               while k3>0;
%                                   k2=k3;
%                                   while k2>0;
%                                       k1=k2;
%                                       while k1>0;
%                                           Sb=Sb+1;
%                                           k1=k1-1;
%                                       end
%                                       k2=k2-1;
%                                   end
%                                   k3=k3-1;
%                               end
%                               k4=k4-1;
%                           end
%                           k5=k5-1;
%                       end
%                       k6=k6-1;
%                   end
%                   k7=k7-1;
%               end
                 Sb=50388;
 
     end
 
     if case_no>=7
         k7=12;
         while k7>(14-initial_state(1,7))
             k6=k7;
             while k6>0
                 k5=k6;
                 while k5>0
                     k4=k5;
                     while k4>0
                         k3=k4;
                         while k3>0
                             k2=k3;
                             while k2>0
                                 k1=k2;
                                 while k1>0
                                     S7=S7+1;
                                     k1=k1-1;
                                 end
                                 k2=k2-1;
                             end
                             k3=k3-1;
                         end
                         k4=k4-1;
                     end
                     k5=k5-1;
                 end
                 k6=k6-1;
             end
             k7=k7-1;
         end
     end
     if case_no>=6
         k6=14-initial_state(1,7)-(initial_state(1,7)==1);
         while k6>(14-initial_state(1,6))
             k5=k6;
             while k5>0
                 k4=k5;
                 while k4>0
                     k3=k4;
                     while k3>0
                         k2=k3;
                         while k2>0
                             k1=k2;
                             while k1>0
                                 S6=S6+1;
                                 k1=k1-1;
                             end
                             k2=k2-1;
                         end
                         k3=k3-1;
                     end
                     k4=k4-1;
                 end
                 k5=k5-1;
             end
             k6=k6-1;
         end
     end
     if case_no>=5
         k5=14-initial_state(1,6)-(initial_state(1,6)==1);
         while k5>(14-initial_state(1,5))
             k4=k5;
             while k4>0
                 k3=k4;
                 while k3>0
                     k2=k3;
                     while k2>0
                         k1=k2;
                         while k1>0
                             S5=S5+1;
                             k1=k1-1;
                         end
                         k2=k2-1;
                     end
                     k3=k3-1;
                 end
                 k4=k4-1;
             end
             k5=k5-1;
         end
     end
     if case_no>=4
         k4=14-initial_state(1,5)-(initial_state(1,5)==1);
         while k4>(14-initial_state(1,4))
             k3=k4;
             while k3>0
                 k2=k3;
                 while k2>0
                     k1=k2;
                     while k1>0
                         S4=S4+1;
                         k1=k1-1;
                     end
                     k2=k2-1;
                 end
                 k3=k3-1;
             end
             k4=k4-1;
         end
     end
     if case_no>=3
         k3=14-initial_state(1,4)-(initial_state(1,4)==1);
         while k3>(14-initial_state(1,3))
             k2=k3;
             while k2>0
                 k1=k2;
                 while k1>0
                     S3=S3+1;
                     k1=k1-1;
                 end
                 k2=k2-1;
             end
             k3=k3-1;
         end
     end
     if case_no>=2
         k2=14-initial_state(1,3)-(initial_state(1,3)==1);
         while k2>(14-initial_state(1,2))
             k1=k2;
             while k1>0
                 S2=S2+1;
                 k1=k1-1;
             end
             k2=k2-1;
         end
     end
     if case_no>=1
         S1=initial_state(1,1)-initial_state(1,2)+(max_flag==0)-(case_no==1);
     end
 
     State_index=Sb+S7+S6+S5+S4+S3+S2+S1;
 
     Transition_matrix(State_index,i)=relative_probabilities(a,1)/sum(relative_probabilities);
 
     a=a+1;
 end   
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 i=i+1;
end

Part a)

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.

...

---

% This script requires Permutation_Generator and Transition_Matrix_Generator to be run first, generating
% the matrix State_list containing the possible states and the transition matrix Transition_matrix. Dimensions of State_list are
% 50388 x 13, with each row containing a distinct state; dimensions of Transition_matrix are 50388 x 50388 (sparse).
Part_a_convert=zeros(13,50388); % Matrix to convert state probabilities to stat probabilities
Part_b_convert=zeros(43,50388); % Matrix to convert state probabilities to net stat boost probabilities
Turn_0=zeros(50388,1); % Vector to store initial state probability
Turn_0(48673,1)=1; % Set probability of all stats at +0 ([7 7 7 7 7 7 7]) at 100% for the initial state
%Turn_0(1,1)=1; % Set probability of all stats at -6 ([1 1 1 1 1 1 1]) at 100% for the initial state
Turn_n=zeros(50388,1); % Vector to store iterated state probabilities
Part_a=zeros(51,13); % output probabilty matrix for stats after n turns
Part_b=zeros(51,45); % output probability matrix for Stored Power after n turns
% Calculate conversion matrix for part a
m=1;
while m<50389
  n=1;
  while n<14
      Part_a_convert(n,m)=(sum(State_list(m,:)==n))/7;
      n=n+1;
  end
  m=m+1;
end
% Calculate output for part a
Part_a(1,:)=(Part_a_convert*Turn_0)';
k=2;
while k<52
  if k==2
      Turn_n=Transition_matrix*Turn_0;
  else
      Turn_n=Transition_matrix*Turn_n;
  end
  Part_a(k,:)=(Part_a_convert*Turn_n)';
  k=k+1;
end

As for results, the full probability distribution is shown in the spoiler below.

...

---

  Turns        -6        -5        -4        -3       -2        -1      +0           +1        +2        +3        +4       +5       +6
0         0         0         0         0         0         0    1.0000         0         0         0         0         0         0
1.0000         0         0         0         0         0    0.1429    0.7143         0    0.1429         0         0         0         0
2.0000         0         0         0         0    0.0204    0.2041    0.5102    0.0408    0.2041         0    0.0204         0         0
3.0000         0         0         0    0.0029    0.0437    0.2187    0.3732    0.0875    0.2187    0.0087    0.0437         0    0.0029
4.0000         0         0    0.0004    0.0083    0.0624    0.2098    0.2851    0.1251    0.2110    0.0250    0.0625    0.0017    0.0087
5.0000         0    0.0001    0.0015    0.0148    0.0744    0.1913    0.2302    0.1498    0.1954    0.0447    0.0751    0.0063    0.0164
6.0000    0.0000    0.0003    0.0032    0.0211    0.0804    0.1711    0.1961    0.1627    0.1793    0.0643    0.0827    0.0136    0.0253
7.0000    0.0000    0.0006    0.0052    0.0264    0.0818    0.1525    0.1745    0.1666    0.1657    0.0814    0.0871    0.0232    0.0350
8.0000    0.0001    0.0012    0.0074    0.0303    0.0803    0.1368    0.1598    0.1648    0.1553    0.0949    0.0897    0.0343    0.0452
9.0000    0.0003    0.0019    0.0095    0.0329    0.0771    0.1240    0.1487    0.1596    0.1478    0.1048    0.0917    0.0459    0.0558
10.0000    0.0005    0.0026    0.0114    0.0343    0.0733    0.1136    0.1396    0.1529    0.1422    0.1115    0.0938    0.0575    0.0667
11.0000    0.0008    0.0034    0.0128    0.0347    0.0693    0.1050    0.1315    0.1458    0.1379    0.1156    0.0962    0.0687    0.0781
12.0000    0.0011    0.0042    0.0139    0.0345    0.0653    0.0978    0.1241    0.1387    0.1342    0.1179    0.0990    0.0793    0.0899
13.0000    0.0014    0.0048    0.0147    0.0339    0.0616    0.0915    0.1171    0.1320    0.1308    0.1189    0.1020    0.0893    0.1020
14.0000    0.0018    0.0054    0.0152    0.0329    0.0581    0.0859    0.1106    0.1258    0.1274    0.1190    0.1050    0.0986    0.1144
15.0000    0.0022    0.0058    0.0154    0.0318    0.0549    0.0807    0.1044    0.1200    0.1240    0.1186    0.1081    0.1073    0.1269
16.0000    0.0025    0.0060    0.0154    0.0305    0.0518    0.0760    0.0987    0.1145    0.1205    0.1178    0.1110    0.1156    0.1397
17.0000    0.0028    0.0062    0.0152    0.0292    0.0489    0.0715    0.0933    0.1093    0.1170    0.1168    0.1137    0.1235    0.1525
18.0000    0.0030    0.0063    0.0149    0.0278    0.0462    0.0674    0.0882    0.1044    0.1135    0.1157    0.1161    0.1311    0.1654
19.0000    0.0032    0.0062    0.0146    0.0264    0.0436    0.0634    0.0834    0.0997    0.1100    0.1144    0.1183    0.1384    0.1783
20.0000    0.0033    0.0061    0.0141    0.0250    0.0411    0.0597    0.0789    0.0952    0.1066    0.1130    0.1203    0.1455    0.1913
21.0000    0.0033    0.0060    0.0136    0.0236    0.0387    0.0561    0.0746    0.0909    0.1031    0.1115    0.1220    0.1523    0.2042
22.0000    0.0033    0.0057    0.0130    0.0223    0.0364    0.0528    0.0705    0.0867    0.0997    0.1099    0.1236    0.1590    0.2171
23.0000    0.0032    0.0055    0.0124    0.0209    0.0341    0.0496    0.0666    0.0826    0.0963    0.1082    0.1249    0.1654    0.2301
24.0000    0.0031    0.0052    0.0117    0.0196    0.0320    0.0465    0.0629    0.0787    0.0930    0.1065    0.1260    0.1717    0.2431
25.0000    0.0030    0.0049    0.0111    0.0183    0.0299    0.0435    0.0593    0.0749    0.0897    0.1046    0.1269    0.1778    0.2560
26.0000    0.0028    0.0046    0.0104    0.0170    0.0279    0.0407    0.0559    0.0712    0.0864    0.1027    0.1277    0.1836    0.2690
27.0000    0.0026    0.0042    0.0096    0.0158    0.0260    0.0380    0.0526    0.0676    0.0832    0.1007    0.1282    0.1894    0.2821
28.0000    0.0024    0.0039    0.0089    0.0146    0.0241    0.0354    0.0494    0.0641    0.0800    0.0986    0.1286    0.1949    0.2952
29.0000    0.0022    0.0035    0.0082    0.0134    0.0223    0.0329    0.0463    0.0607    0.0768    0.0964    0.1288    0.2002    0.3083
30.0000    0.0019    0.0032    0.0075    0.0123    0.0206    0.0305    0.0433    0.0574    0.0736    0.0941    0.1288    0.2054    0.3216
31.0000    0.0017    0.0029    0.0068    0.0112    0.0189    0.0282    0.0404    0.0541    0.0704    0.0917    0.1286    0.2103    0.3349
32.0000    0.0015    0.0025    0.0061    0.0101    0.0173    0.0260    0.0376    0.0509    0.0672    0.0893    0.1282    0.2151    0.3483
33.0000    0.0013    0.0022    0.0054    0.0091    0.0157    0.0238    0.0348    0.0478    0.0641    0.0867    0.1276    0.2197    0.3618
34.0000    0.0011    0.0019    0.0048    0.0081    0.0142    0.0217    0.0322    0.0447    0.0609    0.0841    0.1268    0.2241    0.3754
35.0000    0.0009    0.0016    0.0042    0.0072    0.0128    0.0198    0.0296    0.0417    0.0577    0.0813    0.1258    0.2282    0.3892
36.0000    0.0007    0.0014    0.0036    0.0063    0.0114    0.0179    0.0271    0.0387    0.0545    0.0785    0.1245    0.2322    0.4032
37.0000    0.0006    0.0011    0.0030    0.0055    0.0101    0.0161    0.0247    0.0358    0.0514    0.0755    0.1230    0.2359    0.4173
38.0000    0.0004    0.0009    0.0025    0.0047    0.0088    0.0143    0.0224    0.0330    0.0482    0.0725    0.1213    0.2393    0.4315
39.0000    0.0003    0.0007    0.0020    0.0039    0.0076    0.0126    0.0202    0.0303    0.0451    0.0694    0.1194    0.2425    0.4459
40.0000    0.0002    0.0005    0.0016    0.0032    0.0064    0.0110    0.0180    0.0276    0.0419    0.0662    0.1172    0.2455    0.4605
41.0000    0.0002    0.0004    0.0012    0.0026    0.0054    0.0095    0.0159    0.0249    0.0388    0.0628    0.1148    0.2482    0.4753
42.0000    0.0001    0.0003    0.0009    0.0021    0.0044    0.0080    0.0138    0.0223    0.0356    0.0594    0.1121    0.2506    0.4903
43.0000    0.0001    0.0002    0.0007    0.0016    0.0035    0.0067    0.0118    0.0197    0.0324    0.0558    0.1091    0.2528    0.5056
44.0000    0.0000    0.0001    0.0005    0.0012    0.0028    0.0054    0.0100    0.0171    0.0292    0.0521    0.1058    0.2545    0.5212
45.0000    0.0000    0.0001    0.0004    0.0009    0.0021    0.0043    0.0082    0.0147    0.0259    0.0481    0.1021    0.2558    0.5373
46.0000    0.0000    0.0001    0.0002    0.0006    0.0016    0.0034    0.0067    0.0124    0.0227    0.0440    0.0979    0.2566    0.5538
47.0000    0.0000    0.0000    0.0002    0.0005    0.0012    0.0026    0.0053    0.0102    0.0195    0.0397    0.0931    0.2567    0.5709
48.0000    0.0000    0.0000    0.0001    0.0003    0.0008    0.0019    0.0042    0.0083    0.0165    0.0354    0.0879    0.2561    0.5884
49.0000    0.0000    0.0000    0.0001    0.0002    0.0006    0.0014    0.0032    0.0066    0.0137    0.0311    0.0822    0.2546    0.6063
50.0000    0.0000    0.0000    0.0000    0.0001    0.0004    0.0010    0.0024    0.0051    0.0112    0.0269    0.0760    0.2521    0.6246

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

Note here that they agree well with ganj4lF's simulations. I'd like to thank him for that exact reason, actually - instead of having to ponder endlessly how accurate my results were, his work provided me with a reference point to show that the numbers in my output were probably correct.


Part b)

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.

...

---

% Calculate conversion matrix for part b
m=1;
while m<50389
  n=1;
  while n<44
      Part_b_convert(n,m)=((sum((State_list(m,:)>7).*(State_list(m,:)-7)))==(n-1));
      n=n+1;
  end
  m=m+1;
end
% Calculate output for part b
Part_b(1,1:43)=(Part_b_convert*Turn_0)';
k=2;
while k<52
  if k==2
      Turn_n=Transition_matrix*Turn_0;
  else
      Turn_n=Transition_matrix*Turn_n;
  end
  Part_b(k,1:43)=(Part_b_convert*Turn_n)';
  k=k+1;
end
% Add average boosts and SP BP columns to part b output
k=1;
while k<52
  m=1;
  while m<44
      Part_b(k,44)=Part_b(k,44)+(m-1)*Part_b(k,m);
      m=m+1;
  end
  Part_b(k,45)=20*(1+Part_b(k,44));
  k=k+1;
end

{Having difficulties getting this one to display due to the size of the data set - will fix later}

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

Conclusion

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