Algorithm Fourth Edition: 1.1.27 binomial distribution, from 10 to 10 billion@ TOC
Recursive method to realize binomial distribution operation, O(2^N)
The recursive algorithm is relatively simple. It is very suitable for small parameters.
But the problem is also obvious. N * k > 100 is very slow at the beginning. When > 1000, the program basically crashes. The algorithm efficiency is O(2^N), which can not be used in practical development.
double binomial(int N, int k, double p)
{
if (N == 0 && k == 0)
{
return 1;
}
if (N < 0 || k < 0)
{
return 0;
}
return (1.0 - p) * binomial(N - 1, k, p) + p * binomial(N - 1, k - 1, p);
}
Optimize by saving the calculation results to the array, O(N)
It is also recursive, but saving and reusing the calculated results immediately reduces the complexity and efficiency O(N).
At this time, large number operation can be carried out, for example, n * k = 100000000, which can be used in practice, and it is of little significance.
struct bino
{
double qbinomial(int N, int k, double p)
{
if (N == 0 && k == 0)
{
return 1;
}
if (N < 0 || k < 0)
{
return 0;
}
vd = std::vector<std::vector<double>>((N + 1), std::vector<double>((k + 1), -1));
return qbinomial(0, N, k, p);
}
private:
double qbinomial(int f, int N, int k, double p)
{
if (N == 0 && k == 0)
{
return 1;
}
if (N < 0 || k < 0)
{
return 0;
}
if (vd[N][k] < 0)
{
vd[N][k] = (1.0 - p) * qbinomial(0, N - 1, k, p) + p * qbinomial(0, N - 1, k - 1, p);
}
double result = vd[N][k];
vd.clear();
return result;
}
std::vector<std::vector<double>> vd;
};
Pursue the best and seek a breakthrough
With optimized recursion, the order of magnitude plus one, the program crashes. In order to improve the stability, we abandon the recursion method and use the loop, but we need to have a deep understanding of the algorithm.
Recursion is pushed back and forward, which still takes up too much computation. The stability is improved by pushing directly from the array loop. At the level of N * k =500 million, it does not collapse, but the speed does not increase significantly.
double qqbinomial(int N, int k, double p) //Fast binomial distribution; Building a matrix data structure
{
if (N == 0 && k == 0)
{
return 1;
}
if (N < 0 || k < 0)
{
return 0;
}
std::vector<std::vector<double>> qvd =
std::vector<std::vector<double>>
((N + 1), std::vector<double>((k + 1), 0));
double q = 1;
for (size_t i = 0; i != N + 1; ++i)
{
qvd[i][0] = q;
q *= 1 - p;
}
double pp = 1 - p;
for (size_t i = 1; i != N + 1; ++i)
{
if (i < k) //When I > k, the extra values are 0. Make judgment to avoid repeated calculation
{
for (size_t j = 1; j != i + 1; ++j)
{
qvd[i][j] = pp * qvd[i - 1][j] + p * qvd[i - 1][j - 1];
}
}
else
{
for (size_t j = 1; j != k + 1; ++j)
{
qvd[i][j] = pp * qvd[i - 1][j] + p * qvd[i - 1][j - 1];
}
}
}
double result = qvd[N][k];
qvd.clear();
return result;
}
I want ultimate stability
How can I reach the int limit of C + +? When n * k = 1 billion, there is not enough memory and still crashes, but this has not reached the upper limit of int integer.
Can I reach the upper limit? Probably not, but try.
Reviewing the circular array algorithm, it is found that the number of elements in the array is N * k. when the number exceeds 1 billion, 16g of memory directly crashes. But do I need such a large array?
The problem lies in the data structure. When calculating the value, the algorithm only needs two numbers [n - 1] [k] and [n - 1] [K - 1], which need K numbers. Therefore, the STD:: deque < STD:: vector < double > > data structure can be adopted. Every time a row of array is calculated and inserted, a number group will be removed from the top of the queue. At this time, the occupied memory will not exceed K double memory.
At the same time, according to the algorithm, when k > N, the result is 0, which can improve the efficiency of the algorithm to O(N) or O(1), and the memory space is not greater than the memory occupied by K double types.
At this time, the order of magnitude N * K can reach 10 billion, theoretically breaking through the int upper limit (2 billion), but it is meaningless.
Here are the codes:
double binomial(int N, int k, double p) //Stabilize binomial distribution and construct queue
{
if (N == 0 && k == 0)
{
return 1;
}
if (N < 0 || k < 0)
{
return 0;
}
if (k > N)
{
double result = 0;
return result;
}
std::deque<std::vector<double>> deqVecDou;
deqVecDou.push_back(std::vector<double>{1, 0});
double q = 1;
std::vector<double> lieOne(N + 1);
for (size_t i = 0; i != N + 1; ++i)
{
lieOne[i] = q;
q *= 1 - p;
}
double pp = 1 - p;
for (size_t i = 1; i != N + 1; ++i)
{
if (i < k)
{
deqVecDou.push_back(std::vector<double>(i + 2, 0));
deqVecDou[1][0] = lieOne[i];
for (size_t j = 1; j != i + 1; ++j)
{
deqVecDou[1][j] = pp * deqVecDou[0][j] + p * deqVecDou[0][j - 1];
}
deqVecDou.pop_front();
}
else
{
deqVecDou.push_back(std::vector<double>(k + 2, 0));
deqVecDou[1][0] = lieOne[i];
for (size_t j = 1; j != k + 1; ++j)
{
deqVecDou[1][j] = pp * deqVecDou[0][j] + p * deqVecDou[0][j - 1];
}
deqVecDou.pop_front();
}
}
double result = deqVecDou[0][k];
deqVecDou.clear();
lieOne.clear();
return result;
}
Without data structure, the algorithm is basically a castle in the air.
I spent two days thinking about and improving the small problem in Section 1 of Chapter 1 of the fourth edition of the algorithm. Although it was a blind adjustment, because it had exceeded the actual demand in the end, and the data accuracy of double had been far left behind, it was still meaningful.
Advancing from O(2^N) to O(N) requires only the most basic array. The memory occupation is from N ^ 2 to N. you only need to set an array in the queue, and then delete it at the top of the loop. I deeply realize that the attack of the algorithm is at the rolling level. C + +, which is famous for pursuing efficiency, is basically a joke without the blessing of the algorithm.
For algorithms, see here:
The recursive algorithm of binomial distribution (1.1.27 Binomial distribution) is improved based on array