Join Date: Jun 2012
Posts: 25
Queensland, Australia
|
Problem 1
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.
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.
...
Code:
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.
...
Code:
% 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.
...
Code:
% 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.
However, I will highlight this particular section, the results specifically for turn 50:
Code:
-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.
...
Code:
% 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
Code:
{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.
Code:
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
Last edited by Tsaeb XIII; Jan 1st, 2013 at 11:47:08 PM.
|