[PYTHON] Probability of getting the highest and lowest turnip prices in Atsumori

Overview

The selling price of turnips (weekday prices at Tanuki Shoten) in "Animal Crossing: New Clothes" is said to fluctuate between 9 and 660 bells </ b>, with the highest price of 660 bells and the lowest price of 9 bells. It is said that it is very rare to come out. Therefore, let's calculate the probability that the highest and lowest prices will appear in a certain week.

Premise

  • The probability to be calculated is a value that converges when the number of gameplay weeks approaches infinity.
  • It is assumed that the analysis result of the following turnip value determination algorithm (as of July 13, 2020) is correct. https://gist.github.com/Treeki/85be14d297c80c8b3c0a76375743325b
  • The act of changing the time in the game from the actual one (so-called time operation) shall not be performed.
  • The act of finding regularity in pseudo-random numbers and adjusting them to a specific value (so-called random number adjustment) shall not be performed.

Turnip price fluctuation pattern

The fluctuation pattern of turnip price can be roughly divided into the following four types.

  1. Decrease type </ b>: It decreases monotonically and does not exceed the bid price.
  2. Wave </ b>: Repeat up and down several times to keep the price near the bid price.
  3. Large bounce </ b>: After a monotonous decrease, it starts to increase, reaches a large peak, and decreases.
  4. Small bounce </ b>: After a monotonous decrease, it starts to increase, reaches a small peak, and decreases.

The highest price occurs only in "3. Large bounce", and the lowest price occurs only in "4. Small bounce". Therefore, first try to find the probability of each type.

The probability of becoming each type is determined only by the type of the previous week (Markov chain), as shown in the table below. (Analysis result code lines 135-209)

Last week\this week Reduced type Wave shape Bounce large Bounce small
Reduced type 5% 25% 45% 25%
Wave shape 15% 20% 30% 35%
Bounce large 20% 50% 5% 25%
Bounce small 15% 45% 25% 15%

From this table, try to find the probability (steady distribution) that each type converges when the number of gameplay weeks approaches infinity by the method of the following site. (It is necessary to prove that it actually converges, but it is omitted here.) https://withcation.com/2020/01/23/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96%E3%81%AE%E5%AE%9A%E5%B8%B8%E5%88%86%E5%B8%83%E3%82%92%E8%A7%A3%E3%81%8F/

Assuming that the probabilities of reduced type, wave type, large bounce, and small bounce are $ a, b, c, and d $, respectively, the following equation holds.

\begin{pmatrix}
a & b & c & d
\end{pmatrix}
\begin{pmatrix}
0.05 & 0.25 & 0.45 & 0.25 \\
0.15 & 0.20 & 0.30 & 0.35 \\
0.20 & 0.50 & 0.05 & 0.25 \\
0.15 & 0.45 & 0.25 & 0.15 \\
\end{pmatrix}
=
\begin{pmatrix}
a & b & c & d
\end{pmatrix}

Transfer

\begin{pmatrix}
a & b & c & d
\end{pmatrix}
\begin{pmatrix}
-0.95 & 0.25 & 0.45 & 0.25 \\
0.15 & -0.80 & 0.30 & 0.35 \\
0.20 & 0.50 & -0.95 & 0.25 \\
0.15 & 0.45 & 0.25 & -0.85 \\
\end{pmatrix}
=
\begin{pmatrix}
0 & 0 & 0 & 0
\end{pmatrix}

At first glance, it seems that the simultaneous equations of 4 equations using 4 kinds of characters are expressed in a matrix, but in fact, if any 3 equations are transformed, the remaining 1 equation will be left, so one equation is missing (4x4 matrix). The determinant of is 0, and the inverse matrix cannot be obtained). Therefore, embed $ a + b + c + d = 1 $ in any column (here, the first column).

\begin{pmatrix}
a & b & c & d
\end{pmatrix}
\begin{pmatrix}
1 & 0.25 & 0.45 & 0.25 \\
1 & -0.80 & 0.30 & 0.35 \\
1 & 0.50 & -0.95 & 0.25 \\
1 & 0.45 & 0.25 & -0.85 \\
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0
\end{pmatrix}
\\
\begin{pmatrix}
a & b & c & d
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 0.25 & 0.45 & 0.25 \\
1 & -0.80 & 0.30 & 0.35 \\
1 & 0.50 & -0.95 & 0.25 \\
1 & 0.45 & 0.25 & -0.85 \\
\end{pmatrix}^{-1}
=
\begin{pmatrix}
0.1476 & 0.3463 & 0.2474 & 0.2588
\end{pmatrix}

The probabilities of each type were calculated. The Python code is shown below.

import numpy as np

A = np.array(
    [[0.05, 0.25, 0.45, 0.25],
    [0.15, 0.20, 0.30, 0.35],
    [0.20, 0.50, 0.05, 0.25],
    [0.15, 0.45, 0.25, 0.15]]
    )

A -= np.identity(4)  #Transition
A[:, 0] = 1  #Condition a+b+c+d=Put 1
B = np.array([1, 0, 0, 0])

print(B @ np.linalg.inv(A))  # [0.1476074  0.34627733 0.24736279 0.25875248]

Probability of reaching the highest value (660 bells)

In order to reach the maximum value (660 bells), in addition to "3. Large bounce" (probability 0.2474), the following conditions are required.

A: The island's bid price (the price you can buy from Uri on Sunday morning) is 90-110 bells, which is 110 bells. </ b> B: At peak times, it is 2 to 6 times the bid price, but this is 6 times. </ b>

For A, it is random (randint on line 123 of the code) and the probability is $ \ frac {1} {21} $. For B, first generate a float type numerical value of 2.0 to 6.0 as a magnification, but the peak price is not exactly 6.0 times because it is the result of (buying price) × (magnification) rounded up to the nearest whole number. There is a possibility that it may be. The magnification is determined based on a random number of type uint32_t (takes an integer value from $ 0 $ to $ 2 ^ {32} $), and the larger the random number value, the larger the magnification. Therefore, the probability of satisfying condition B is obtained by examining the border of uint32_t type random numbers whose peak price is 660 bells using a binary search.

#include <bits/stdc++.h>
using namespace std;

// uint32_A function that returns a float type number based on a t type random number. The larger the random number, the closer to b.
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

//Function to round up after the decimal point
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

int main() {
  uint32_t ll = 0;  //uint32_Lower limit of t
  uint32_t rr = pow(2, 32) - 1;  //uint32_Upper limit of t
  int32_t basePrice = 110;  //Bid price
  uint64_t l, r, m;
  int32_t peakPrice;

  l = 0;
  r = pow(2, 32) - 1;
  m = (l + r) / 2;

  //Binary search
  while (r - l > 1) {
    peakPrice = intceil(randfloat(2.0, 6.0, m) * basePrice);  //Peak price
    if (peakPrice < 660) {l = m;}
    else {r = m;}
    m = (l + r) / 2;
  }

  cout << 0 << "From" << rr << "The random number that takes the value of" << m << "If it is larger, it will be the highest value." << endl;
  cout << "The probability is" << (double)(rr - m) / (double)(rr) << "Is." << endl;
}

output

If the random number that takes a value between 0 and 4294967295 is greater than 4285206015, it is the highest value.
The probability is 0.It is 00227273.

Multiplying all the probabilities found gives $ 0.2474 \ times \ frac {1} {21} \ times 0.002272 = 2.68 \ times 10 ^ {-5} $. Although it is a low probability of 0.00268% </ b>, the number of games sold is <a href="https://www.nintendo.co.jp/ir/pdf/2020/" as of May 7, 2020. 200507_3.pdf "> 11.77 million bottles , so even if 1% of them check the turnip price, about 3 people around the world check the maximum value of 660 bells every week. It becomes the calculation that is being done.

Next, let's calculate the probability of the lowest price.

Probability of the lowest price (9 bells)

In order to reach the lowest price (9 bells), in addition to "4. Bounce small size" (probability 0.2588), the following conditions are required.

A: The island's bid price (the price you can buy from Uri on Sunday morning) is 90-110 bells, which is 90 bells. </ b> B: Of the eight "time to start increasing for the first time" {Monday morning, Monday afternoon, ..., Thursday afternoon}, Monday morning or Thursday afternoon. </ b>

  • In B, the probabilities for Monday morning and Thursday afternoon are the same, and the programs are almost the same, so the conditions for C and D below will be explained for Thursday afternoon. </ b> C: Under these conditions, the turnip price on Monday morning could be 36-81 bells, which is 36 bells. </ b> D: Under these conditions, the turnip price on Thursday afternoon could be 9-65 bells, which is 9 bells. </ b>

Since A and B are selected by random numbers of integers, the probabilities are $ \ frac {1} {21} and \ frac {1} {4} $, respectively. For C, first generate a float type numerical value of 0.4 to 0.9 as a magnification, and multiply it by 90 bells of the bid price to obtain the turnip price on Monday morning. Since the number after the decimal point is rounded up, the magnification must be 0.4 (due to the error, there is no problem if it is slightly higher) in order to reach 36 bells. As mentioned above, the magnification is determined based on the random number of uint32_t type (takes an integer value from $ 0 $ to $ 2 ^ {32} $), but here, the smaller the random number value, the larger the magnification. (Code line 317). After all, it seems good to find the probability of satisfying condition C by examining the border of a random number with a turnip value of 36 bells using a binary search.

#include <bits/stdc++.h>
using namespace std;

// uint32_A function that returns a float type number based on a t type random number. The larger the random number, the closer to b.
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

//Function to round up after the decimal point
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

int main() {
  uint32_t ll = 0;  //uint32_Lower limit of t
  uint32_t rr = pow(2, 32) - 1;  //uint32_Upper limit of t
  int32_t basePrice = 90;  //Bid price
  uint64_t l, r, m;
  int32_t monAmPrice;

  l = 0;
  r = pow(2, 32) - 1;
  m = (l + r) / 2;
  
  //Binary search
  while (r - l > 1) {
    monAmPrice = intceil(randfloat(0.9, 0.4, m) * basePrice);  //Monday morning price
    if (monAmPrice > 36) {l = m;}
    else {r = m;}
    m = (l + r) / 2;
  }

  cout << ll << "From" << rr << "The random number that takes the value of" << m << "If greater than, condition C is satisfied." << endl;
  cout << "The probability is" << (double)(rr - m) / (double)(rr) << "Is." << endl;
}

output

Condition C is satisfied if the random number from 0 to 4294967295 is greater than 4294966783.
The probability is 1.19209e-07.

Due to the nature of rounding up to the nearest whole number, the probability is much lower this time than condition B at the highest value.

Finally, consider the probability of satisfying the condition "D: Thursday afternoon turnip value can be 9-65 bells, which is 9 bells." The multiplier, which is 0.4 as of Monday morning, decreases monotonically over the six periods from Monday afternoon to Thursday afternoon. At this time, in each period, the amount of decrease in the magnification is determined from 0.03 to 0.05 and updated. In order for 90 bells to become 9 bells, the magnification must be $ 0.1 = 0.4-0.05 \ times 6 $, and we will consider the case where the rate of decline is almost the maximum (0.05) for 6 consecutive periods. However, when rounding up the turnip value after the decimal point, the process of "adding 0.9999 and then rounding down the decimals" is performed, so for example, 9.001 bells are rounded up to 10 bells, but 9.0000001 bells become 9 bells. To do. Therefore, it is possible that the rate of decline does not have to be the maximum value for all six times. Even so, under condition C, the threshold value of the random number was obtained by performing a binary search, but that was possible because the random number was generated only once, and this time 6 random numbers are involved, so simply Can not. Therefore, the probability is calculated by the following method.

D-1. Let $ A $ be the amount of decrease in magnification (closest to 0.03 (Note)) when the random number is the highest value $ 2 ^ {32} -1 $. Under the condition that the amount of decline of 5 out of 6 times is $ A $, the threshold value of how many or more of the remaining 1 decrease will be 9 bells is calculated. Let the threshold be $ B $. D-2. For the amount of decrease from $ B $ to $ A $, find the range of the random number to be the amount of decrease. D-3. Search through each of the 6 drops between $ B $ and $ A $, and for those with 9 bells, the number of such cases (multiply the width calculated by D-2 by 6). Add together). D-4. The value divided by $ (2 ^ {32}) ^ 6 $ is the probability of satisfying condition D.

(Note) The reason why it is 0.03 instead of 0.05 is that the analysis uses the operation of "subtract 0.02 and then subtract 0 to 0.03", and it follows that.

Regarding the amount of decline, it is easier to handle by looking at the mantissa part of the float type than displaying the decimal point in decimal numbers, so it is written as such. I referred to the following article to get the hexadecimal representation of the mantissa. https://qiita.com/nia_tn1012/items/d26f0fc993895a09b30b In addition, the value of random numbers, which was previously expressed in decimal, will now be expressed in hexadecimal.

First, perform D-1 by binary search.

#include <bits/stdc++.h>
using namespace std;

// uint32_A function that returns a float type number based on a t type random number. The larger the random number, the closer to b.
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

//Function to round up after the decimal point
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

//A function that returns the mantissa of float type
int printb(float ff) {
    union { float f; int i; } a;
    int i;
    a.f = ff;

    return a.i & 0x7FFFFF;
  }

int main() {
  uint32_t ll = 0;  //uint32_Lower limit of t
  uint32_t rr = pow(2, 32) - 1;  //uint32_Upper limit of t
  int32_t basePrice = 90;  //Bid price
  uint64_t l, r, m;

  float declineRate;
  int32_t thuPmPrice;
  l = 0;
  r = pow(2, 32) - 1;
  m = (l + r) / 2;

  //Binary search
  while (r - l > 1) {
    declineRate = (0.4 - 0.03 * 6 - randfloat(0, 0.02, rr) * 5 - randfloat(0, 0.02, m));
    thuPmPrice = intceil(declineRate * basePrice);  //Thursday afternoon price
    if (thuPmPrice > 9) {l = m;}
    else {r = m;}
    m = (l + r) / 2;
  }
  cout << hex << "The mantissa part of the fall width B is" << printb(randfloat(0, 0.02, m+1)) \
  << "And the value of the smallest random number in the range of decline is" << m+1 << "Is." << endl;
  cout << hex << "The mantissa part of the fall width A is" << printb(randfloat(0, 0.02, rr)) << "Is." << endl;
}

output

The mantissa part of the fall width B is 23d6db, and the smallest random number value in the fall width is ffffb600.
The mantissa part of the fall width A is 23d709.

Therefore, you can search for 23d6db to 23d709, but just in case, check 48 ways of 23d6da to 23d709. (Depending on the calculation order, the mantissa may be 9 bells even if the drop is 23d6da due to rounding error.)

Next, D-2 to D-4 were performed at once. The calculation time took about 10 minutes.

#include <bits/stdc++.h>
using namespace std;

// uint32_A function that returns a float type number based on a t type random number. The larger the random number, the closer to b.
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

//Function to round up after the decimal point
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

//A function that returns the mantissa of float type
int printb(float ff) {
    union { float f; int i; } a;
    int i;
    a.f = ff;
    return a.i & 0x7FFFFF;
  }

//Recursive function that searches all the ways in which the magnification drops
int64_t fullSearch(float declineRateNow, int step, int basePrice, vector<float> declineRates, int c) {
  int64_t ans = 0;
  int32_t thuPmPrice;
  if (step == 6) {
    thuPmPrice = intceil(declineRateNow * basePrice); //Check if the price on Thursday afternoon is 9 bells
    if (thuPmPrice == 9) {return 1;}
    else {return 0;}
  }
  else {
    declineRateNow -= 0.03;
    for (int i=0; i<c; i++) {
      ans += fullSearch(declineRateNow - declineRates.at(i), step + 1, basePrice, declineRates, c);
    }
  }
  return ans;
}

int main() {
  uint32_t m = pow(2, 32) - 1;
  float declineRate = randfloat(0, 0.02, m);
  float declineRate_b = randfloat(0, 0.02, m);
  uint32_t l, r;
  int c = 0;
  vector<float> declineRates(50);
  
  // D-2
  r = m;
  cout << "No. |Random number lower limit|Random number upper limit|width|Mantissa part of the amount of decline" << endl;
  while (true) {
    declineRate_b = declineRate;
    declineRate = randfloat(0, 0.02, m);
    if (printb(declineRate) < printb(declineRate_b)) {
      l = m + 1;
      cout << hex << setfill('0') << setw(2) << c << " | " << l << " | " << r << " | " \
      << r-l+1 << " | " << printb(declineRate_b) << endl;
      declineRates.at(c) = declineRate_b;
      r = m;
      c++;
      if (printb(declineRate) < 0x23d6da) {break;}
    }
    m--;
  }

  // D-3
  int32_t basePrice = 90;
  int64_t count9bell = fullSearch(0.4, 0, basePrice, declineRates, c);

  // D-4
  // (2^9)^After multiplying by 6(2^32)^23 by 2 instead of dividing by 6*Divided 6 times
  double prob_d = count9bell;
  for (int i=0; i<23*6; i++) {
    prob_d /= 2;
  }
  cout << dec << "37^Of the 6 ways to lower the magnification, the turnip price is 9 bells." << count9bell << "It's a street." << endl;
  cout << "In that case, condition D is met and the probability is" << prob_d << "Is." << endl;
}

First, as D-2, we calculated the range of the fall in the 48 ways of the mantissa part from 23d6da to 23d709, and the range of the fall. The following output was obtained.

No. |Random number lower limit|Random number upper limit|width|Mantissa part of the amount of decline
00 | fffffe00 | ffffffff | 200 | 23d709
01 | fffffc00 | fffffdff | 200 | 23d707
02 | fffffa00 | fffffbff | 200 | 23d706
03 | fffff800 | fffff9ff | 200 | 23d705
04 | fffff600 | fffff7ff | 200 | 23d704
05 | fffff400 | fffff5ff | 200 | 23d702
06 | fffff200 | fffff3ff | 200 | 23d701
07 | fffff000 | fffff1ff | 200 | 23d700
08 | ffffee00 | ffffefff | 200 | 23d6fe
09 | ffffec00 | ffffedff | 200 | 23d6fd
0a | ffffea00 | ffffebff | 200 | 23d6fc
0b | ffffe800 | ffffe9ff | 200 | 23d6fb
0c | ffffe600 | ffffe7ff | 200 | 23d6f9
0d | ffffe400 | ffffe5ff | 200 | 23d6f8
0e | ffffe200 | ffffe3ff | 200 | 23d6f7
0f | ffffe000 | ffffe1ff | 200 | 23d6f6
10 | ffffde00 | ffffdfff | 200 | 23d6f4
11 | ffffdc00 | ffffddff | 200 | 23d6f3
12 | ffffda00 | ffffdbff | 200 | 23d6f2
13 | ffffd800 | ffffd9ff | 200 | 23d6f0
14 | ffffd600 | ffffd7ff | 200 | 23d6ef
15 | ffffd400 | ffffd5ff | 200 | 23d6ee
16 | ffffd200 | ffffd3ff | 200 | 23d6ed
17 | ffffd000 | ffffd1ff | 200 | 23d6eb
18 | ffffce00 | ffffcfff | 200 | 23d6ea
19 | ffffcc00 | ffffcdff | 200 | 23d6e9
1a | ffffca00 | ffffcbff | 200 | 23d6e7
1b | ffffc800 | ffffc9ff | 200 | 23d6e6
1c | ffffc600 | ffffc7ff | 200 | 23d6e5
1d | ffffc400 | ffffc5ff | 200 | 23d6e4
1e | ffffc200 | ffffc3ff | 200 | 23d6e2
1f | ffffc000 | ffffc1ff | 200 | 23d6e1
20 | ffffbe00 | ffffbfff | 200 | 23d6e0
21 | ffffbc00 | ffffbdff | 200 | 23d6de
22 | ffffba00 | ffffbbff | 200 | 23d6dd
23 | ffffb800 | ffffb9ff | 200 | 23d6dc
24 | ffffb600 | ffffb7ff | 200 | 23d6db

All numbers in the table are hexadecimal. For example, looking at the last line, if (random numbers taking values from 00000000 to ffffffff) take 0x200 values from ffffb600 to ffffb7ff, the mantissa part of the fall will be 23d6db. From this table, we can see the following.

  • No matter what the random number is, the amount of decline will not be the mantissa part such as 23d708, 23d703, .... (Although we have examined 48 ways, the number of rows in the table is 0x25, that is, 37 rows, so we can see that there are 11 such declines.)
  • For possible drop widths, the width of the random number range is all 0x200 or $ 2 ^ 9 = 512 $. (Actually, this can be seen from the code of the randfloat function)

Therefore, for D-3 and D-4, we searched all 6 times of decline in these 37 ways, and counted the number of times that the turnip value became 9 bells in the $ 37 ^ 6 $ pattern of decline, and the number of times Multiplied by $ (2 ^ 9) ^ 6 $ and then divided by $ (2 ^ {32}) ^ 6 $, but the condition "D: Thursday afternoon turnip price may be 9-65 bells" However, this is the probability of satisfying "9 bells." As a result of calculating this, the following output was obtained.

37^Of the six "how to reduce the magnification", 8689156 have a turnip price of 9 bells.
In that case, condition D is met and the probability is 2.49367e-It is 035.

Multiplying the probability of the fluctuation pattern (small bounce) and the probabilities of conditions A to D, $ 0.2588 \ times \ frac {1} {21} \ times \ frac {1} {4} \ times (1.192 \ times 10 ^ {-7}) \ times (2.494 \ times 10 ^ {-35}) = 9.16 \ times 10 ^ {-45} $. This is an astronomical low probability. Even if all 7.7 billion people living on Earth have this game and check the turnip price every week from the birth of the universe 13.8 billion years ago to the present, even once they see "1 turnip 9 bells". The probability of being done is $ 5.08 \ times 10 ^ {-23} $. </ b> It should be noted that this is the probability when the source code of the analysis result of the turnip value determination algorithm described in "Premise" is completely the same as that inside the game. If the value added by the intceil function that rounds up after the decimal point is 0.999999 instead of 0.99999, the probability of the lowest price will be lower, and even if the float type is double type, the probability will change.

Summary

  • The probability that the turnip price will reach the highest value of 660 bells is $ 2.68 \ times 10 ^ {-5} $, and it can be said that it appears somewhere in the world every week.
  • The probability that the turnip price will be the lowest 9 bells is $ 9.16 \ times 10 ^ {-45} $, which is unlikely to happen.

reference

Animal Crossing: New Year's Forest Turnip Price Prediction Tool | hyperWiki When you enter the buy price or sell price of the turnip, the predicted value of the turnip price and the probability that the turnip price is each fluctuation pattern are displayed.

AtCoder code test Code test page of the competitive programming site AtCoder. Even if Japanese is used for the standard output of C ++, it is output without garbled characters.

Change log

  • 2020/07/14 Corrected the probability of the lowest price, corrected typographical errors (wrong: rate of decline → correct: width of decline)

Recommended Posts