本部落格所有文章分類的總目錄:【總目錄】本部落格博文總目錄-實時更新
開源Math.NET基礎數學類庫使用總目錄:【目錄】開源Math.NET基礎數學類庫使用總目錄
前言
真正意義上的隨機數(或者隨機事件)在某次產生過程中是按照實驗過程中表現的分佈概率隨機產生的,其結果是不可預測的,是不可見的。而計算機中的隨機函式是按照一定演算法模擬產生的,其結果是確定的,是可見的。我們可以這樣認為這個可預見的結果其出現的概率是100%。所以用計算機隨機函式所產生的“隨機數”並不隨機,是偽隨機數。偽隨機數的作用在開發中的使用非常常見,因此.NET在System名稱空間,提供了一個簡單的Random隨機數生成型別。但這個型別並不能滿足所有的需求,本節開始就將陸續介紹Math.NET中有關隨機數的擴充套件以及其他偽隨機生成演算法編寫的隨機數生成器。
今天要介紹的是Math.NET中擴充套件的其他隨機數生成演算法。
如果本文資源或者顯示有問題,請參考 本文原文地址:http://www.cnblogs.com/asxinyu/p/4301555.html
http://zh.wikipedia.org/wiki/隨機數
隨機數是專門的隨機試驗的結果。
在統計學的不同技術中需要使用隨機數,比如在從統計總體中抽取有代表性的樣本的時候,或者在將實驗動物分配到不同的試驗組的過程中,或者在進行蒙特卡羅模擬法計算的時候等等。產生隨機數有多種不同的方法。這些方法被稱為隨機數生成器。隨機數最重要的特性是它在產生時後面的那個數與前面的那個數毫無關係。
真正的隨機數是使用物理現象產生的:比如擲錢幣、骰子、轉輪、使用電子元件的噪音、核裂變等等。這樣的隨機數生成器叫做物理性隨機數生成器,它們的缺點是技術要求比較高。在實際應用中往往使用偽隨機數就足夠了。這些數列是“似乎”隨機的數,實際上它們是通過一個固定的、可以重複的計算方法產生的。它們不真正地隨機,因為它們實際上是可以計算出來的,但是它們具有類似於隨機數的統計特徵。這樣的生成器叫做偽隨機數生成器。在真正關鍵性的應用中,比如在密碼學中,人們一般使用真正的隨機數。
1.Math.NET的其他隨機數生成演算法
Math.NET在MathNet.Numerics.Random名稱空間下包括了若干個其他隨機數生成演算法,基本介紹如下:
1.Mcg31m1類 與 Mcg59 類,都屬於 矩陣同餘發生器,矩陣同餘發生器是乘同餘線性發生器的一個推廣;2者的引數有些差別;
3.MersenneTwister,Mersenne Twister演算法譯為馬特賽特旋轉演演算法,是偽隨機數發生器之一,其主要作用是生成偽隨機數。此演算法是Makoto Matsumoto (松本)和Takuji Nishimura (西村)於1997年開發的,基於有限二進位制欄位上的矩陣線性再生。可以快速產生高質量的偽隨機數,修正了古老隨機數產生演算法的很多缺陷。 Mersenne Twister這個名字來自週期長度通常取Mersenne質數這樣一個事實。常見的有兩個變種Mersenne Twister MT19937和Mersenne Twister MT19937-64。Math.NET實現的版本是前者(Mersenne Twister MT19937)。介紹
4.Mrg32k3a,又叫 素數模乘同餘法,素數模乘同餘發生器是提出的一個統計性質較好和週期較大的發生器,目前是使用最廣的一種均勻隨機數發生器下面給出兩組經檢驗統計性質是良好的素數模乘同餘發生器;
5.Palf,是一個並行加法滯後的斐波那契偽隨機數字生成器。
6.WH1982與WH2006類,是乘線性同餘法,也叫積式發生器,Math.NET實現的2個版本方便是1982和2006,代表作者釋出該演算法論文的2個時間,可以參考作者的2篇論文:
1.Wichmann, B. A. & Hill, I. D. (1982), "Algorithm AS 183:An efficient and portable pseudo-random number generator". Applied Statistics 31 (1982) 188-190
2. Wichmann, B. A. & Hill, I. D. (2006), "Generating good pseudo-random numbers".Computational Statistics & Data Analysis 51:3 (2006) 1614-1622
8.Xorshift類,實現的是George Marsaglia,在論文“Xorshift RNGs”提出的一個演算法,原理可以參考這篇論文:http://www.jstatsoft.org/v08/i14/paper;
2.Math.NET擴充套件隨機數生成演算法的實現
上面已經已經對Math.NET擴充套件的幾個隨機數生成演算法進行了介紹,對於一般人來說,直接用就可以了,但對於特殊的人來說,可能要用到其中一種,可以直接使用C#進行呼叫即可,當然為了節省大家的時間,這裡對Math.NET的實現做一個簡單的介紹,這樣大家可以更加快速的擴充套件自己的演算法,同時也可以瞭解實現的原理,對Math.NET有一個更加深入的瞭解。
Math.NET對隨機數的擴充套件也是借用了System.Random類,在它的基礎上實現了一個隨機數發生器的基類:RandomSource,所有的擴充套件演算法都實現該類,這樣使用RandomSource就非常方便。RandomSource的結構很簡單,對System.Random進行了簡單的封裝,增加了幾個直接生成其他型別隨機數的方法,並都可以在繼承中使用。其原始碼如下:
1 public abstract class RandomSource : System.Random 2 { 3 readonly bool _threadSafe; 4 readonly object _lock = new object(); 5 6 /// <summary> 7 /// Initializes a new instance of the <see cref="RandomSource"/> class using 8 /// the value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to set whether 9 /// the instance is thread safe or not. 10 /// </summary> 11 protected RandomSource() : base(RandomSeed.Robust()) 12 { 13 _threadSafe = Control.ThreadSafeRandomNumberGenerators; 14 } 15 16 /// <summary> 17 /// Initializes a new instance of the <see cref="RandomSource"/> class. 18 /// </summary> 19 /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param> 20 /// <remarks>Thread safe instances are two and half times slower than non-thread 21 /// safe classes.</remarks> 22 protected RandomSource(bool threadSafe) : base(RandomSeed.Robust()) 23 { 24 _threadSafe = threadSafe; 25 } 26 27 /// <summary> 28 /// Fills an array with uniform random numbers greater than or equal to 0.0 and less than 1.0. 29 /// </summary> 30 /// <param name="values">The array to fill with random values.</param> 31 public void NextDoubles(double[] values) 32 { 33 if (_threadSafe) 34 { 35 lock (_lock) 36 { 37 for (var i = 0; i < values.Length; i++) 38 { 39 values[i] = DoSample(); 40 } 41 } 42 } 43 else 44 { 45 for (var i = 0; i < values.Length; i++) 46 { 47 values[i] = DoSample(); 48 } 49 } 50 } 51 52 /// <summary> 53 /// Returns an infinite sequence of uniform random numbers greater than or equal to 0.0 and less than 1.0. 54 /// </summary> 55 public IEnumerable<double> NextDoubleSequence() 56 { 57 for (int i = 0; i < 64; i++) 58 { 59 yield return NextDouble(); 60 } 61 62 var buffer = new double[64]; 63 while (true) 64 { 65 NextDoubles(buffer); 66 for (int i = 0; i < buffer.Length; i++) 67 { 68 yield return buffer[i]; 69 } 70 } 71 } 72 73 /// <summary> 74 /// Returns a nonnegative random number. 75 /// </summary> 76 /// <returns> 77 /// A 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>. 78 /// </returns> 79 public override sealed int Next() 80 { 81 if (_threadSafe) 82 { 83 lock (_lock) 84 { 85 return (int)(DoSample()*int.MaxValue); 86 } 87 } 88 89 return (int)(DoSample()*int.MaxValue); 90 } 91 92 /// <summary> 93 /// Returns a random number less then a specified maximum. 94 /// </summary> 95 /// <param name="maxValue">The exclusive upper bound of the random number returned.</param> 96 /// <returns>A 32-bit signed integer less than <paramref name="maxValue"/>.</returns> 97 /// <exception cref="T:System.ArgumentOutOfRangeException"><paramref name="maxValue"/> is negative. </exception> 98 public override sealed int Next(int maxValue) 99 { 100 if (maxValue <= 0) 101 { 102 throw new ArgumentException(Resources.ArgumentMustBePositive); 103 } 104 105 if (_threadSafe) 106 { 107 lock (_lock) 108 { 109 return (int)(DoSample()*maxValue); 110 } 111 } 112 113 return (int)(DoSample()*maxValue); 114 } 115 116 /// <summary> 117 /// Returns a random number within a specified range. 118 /// </summary> 119 /// <param name="minValue">The inclusive lower bound of the random number returned.</param> 120 /// <param name="maxValue">The exclusive upper bound of the random number returned. <paramref name="maxValue"/> must be greater than or equal to <paramref name="minValue"/>.</param> 121 /// <returns> 122 /// A 32-bit signed integer greater than or equal to <paramref name="minValue"/> and less than <paramref name="maxValue"/>; that is, the range of return values includes <paramref name="minValue"/> but not <paramref name="maxValue"/>. If <paramref name="minValue"/> equals <paramref name="maxValue"/>, <paramref name="minValue"/> is returned. 123 /// </returns> 124 /// <exception cref="T:System.ArgumentOutOfRangeException"><paramref name="minValue"/> is greater than <paramref name="maxValue"/>. </exception> 125 public override sealed int Next(int minValue, int maxValue) 126 { 127 if (minValue > maxValue) 128 { 129 throw new ArgumentException(Resources.ArgumentMinValueGreaterThanMaxValue); 130 } 131 132 if (_threadSafe) 133 { 134 lock (_lock) 135 { 136 return (int)(DoSample()*(maxValue - minValue)) + minValue; 137 } 138 } 139 140 return (int)(DoSample()*(maxValue - minValue)) + minValue; 141 } 142 143 /// <summary> 144 /// Fills an array with random numbers within a specified range. 145 /// </summary> 146 /// <param name="values">The array to fill with random values.</param> 147 /// <param name="minValue">The inclusive lower bound of the random number returned.</param> 148 /// <param name="maxValue">The exclusive upper bound of the random number returned. <paramref name="maxValue"/> must be greater than or equal to <paramref name="minValue"/>.</param> 149 public void NextInt32s(int[] values, int minValue, int maxValue) 150 { 151 if (_threadSafe) 152 { 153 lock (_lock) 154 { 155 for (var i = 0; i < values.Length; i++) 156 { 157 values[i] = (int)(DoSample()*(maxValue - minValue)) + minValue; 158 } 159 } 160 } 161 else 162 { 163 for (var i = 0; i < values.Length; i++) 164 { 165 values[i] = (int)(DoSample()*(maxValue - minValue)) + minValue; 166 } 167 } 168 } 169 170 /// <summary> 171 /// Returns an infinite sequence of random numbers within a specified range. 172 /// </summary> 173 /// <param name="minValue">The inclusive lower bound of the random number returned.</param> 174 /// <param name="maxValue">The exclusive upper bound of the random number returned. <paramref name="maxValue"/> must be greater than or equal to <paramref name="minValue"/>.</param> 175 public IEnumerable<int> NextInt32Sequence(int minValue, int maxValue) 176 { 177 for (int i = 0; i < 64; i++) 178 { 179 yield return Next(minValue, maxValue); 180 } 181 182 var buffer = new int[64]; 183 while (true) 184 { 185 NextInt32s(buffer, minValue, maxValue); 186 for (int i = 0; i < buffer.Length; i++) 187 { 188 yield return buffer[i]; 189 } 190 } 191 } 192 193 /// <summary> 194 /// Fills the elements of a specified array of bytes with random numbers. 195 /// </summary> 196 /// <param name="buffer">An array of bytes to contain random numbers.</param> 197 /// <exception cref="T:System.ArgumentNullException"><paramref name="buffer"/> is null. </exception> 198 public override void NextBytes(byte[] buffer) 199 { 200 if (buffer == null) 201 { 202 throw new ArgumentNullException("buffer"); 203 } 204 205 if (_threadSafe) 206 { 207 lock (_lock) 208 { 209 for (var i = 0; i < buffer.Length; i++) 210 { 211 buffer[i] = (byte)(((int)(DoSample()*int.MaxValue))%256); 212 } 213 } 214 215 return; 216 } 217 218 for (var i = 0; i < buffer.Length; i++) 219 { 220 buffer[i] = (byte)(((int)(DoSample()*int.MaxValue))%256); 221 } 222 } 223 224 /// <summary> 225 /// Returns a random number between 0.0 and 1.0. 226 /// </summary> 227 /// <returns>A double-precision floating point number greater than or equal to 0.0, and less than 1.0.</returns> 228 protected override sealed double Sample() 229 { 230 if (_threadSafe) 231 { 232 lock (_lock) 233 { 234 return DoSample(); 235 } 236 } 237 238 return DoSample(); 239 } 240 241 /// <summary> 242 /// Returns a random number between 0.0 and 1.0. 243 /// </summary> 244 /// <returns> 245 /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. 246 /// </returns> 247 protected abstract double DoSample(); 248 }
在擴充套件隨機數生成演算法的時候,直接繼承該類即可,看一下Mcg59的實現原始碼:
1 public class Mcg59 : RandomSource 2 { 3 const ulong Modulus = 576460752303423488; 4 const ulong Multiplier = 302875106592253; 5 const double Reciprocal = 1.0/Modulus; 6 ulong _xn; 7 8 /// <summary> 9 /// Initializes a new instance of the <see cref="Mcg59"/> class using 10 /// a seed based on time and unique GUIDs. 11 /// </summary> 12 public Mcg59() : this(RandomSeed.Robust()) 13 { 14 } 15 16 /// <summary> 17 /// Initializes a new instance of the <see cref="Mcg59"/> class using 18 /// a seed based on time and unique GUIDs. 19 /// </summary> 20 /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param> 21 public Mcg59(bool threadSafe) : this(RandomSeed.Robust(), threadSafe) 22 { 23 } 24 25 /// <summary> 26 /// Initializes a new instance of the <see cref="Mcg59"/> class. 27 /// </summary> 28 /// <param name="seed">The seed value.</param> 29 /// <remarks>If the seed value is zero, it is set to one. Uses the 30 /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to 31 /// set whether the instance is thread safe.</remarks> 32 public Mcg59(int seed) 33 { 34 if (seed == 0) 35 { 36 seed = 1; 37 } 38 39 _xn = (uint)seed%Modulus; 40 } 41 42 /// <summary> 43 /// Initializes a new instance of the <see cref="Mcg59"/> class. 44 /// </summary> 45 /// <param name="seed">The seed value.</param> 46 /// <remarks>The seed is set to 1, if the zero is used as the seed.</remarks> 47 /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param> 48 public Mcg59(int seed, bool threadSafe) : base(threadSafe) 49 { 50 if (seed == 0) 51 { 52 seed = 1; 53 } 54 55 _xn = (uint)seed%Modulus; 56 } 57 58 /// <summary> 59 /// Returns a random number between 0.0 and 1.0. 60 /// </summary> 61 /// <returns> 62 /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. 63 /// </returns> 64 protected override sealed double DoSample() 65 { 66 double ret = _xn*Reciprocal; 67 _xn = (_xn*Multiplier)%Modulus; 68 return ret; 69 } 70 71 /// <summary> 72 /// Fills an array with random numbers greater than or equal to 0.0 and less than 1.0. 73 /// </summary> 74 /// <remarks>Supports being called in parallel from multiple threads.</remarks> 75 public static void Doubles(double[] values, int seed) 76 { 77 if (seed == 0) 78 { 79 seed = 1; 80 } 81 82 ulong xn = (uint)seed%Modulus; 83 84 for (int i = 0; i < values.Length; i++) 85 { 86 values[i] = xn*Reciprocal; 87 xn = (xn*Multiplier)%Modulus; 88 } 89 } 90 91 /// <summary> 92 /// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0. 93 /// </summary> 94 /// <remarks>Supports being called in parallel from multiple threads.</remarks> 95 [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] 96 public static double[] Doubles(int length, int seed) 97 { 98 var data = new double[length]; 99 Doubles(data, seed); 100 return data; 101 } 102 103 /// <summary> 104 /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. 105 /// </summary> 106 /// <remarks>Supports being called in parallel from multiple threads, but the result must be enumerated from a single thread each.</remarks> 107 public static IEnumerable<double> DoubleSequence(int seed) 108 { 109 if (seed == 0) 110 { 111 seed = 1; 112 } 113 114 ulong xn = (uint)seed%Modulus; 115 116 while (true) 117 { 118 yield return xn*Reciprocal; 119 xn = (xn*Multiplier)%Modulus; 120 } 121 } 122 }
隨機數的使用大家都很熟練,和Random類差不多,上述擴充套件的演算法中,也包括了很多 靜態方法,可以直接使用。這裡不再舉例說明。
3.資源
原始碼下載:http://www.cnblogs.com/asxinyu/p/4264638.html
如果本文資源或者顯示有問題,請參考 本文原文地址:http://www.cnblogs.com/asxinyu/p/4301555.html