Mersenne Twister (梅森旋转算法）的MATLAB实现

这里写自定义目录标题

• Mersenne Twister (梅森旋转算法）的MATLAB实现
• • Mersenne Twister 简介
  • Mersenne Twister 的使用
  • Mersenne Twister 相关代码

Mersenne Twister (梅森旋转算法）的MATLAB实现

想自己编写一个MATLAB的梅森旋转算法产生随机数的程序，但是找遍网络也没找到拿MATLAB写的，所以试着写一个~

适合学过相关统计知识的朋友

Mersenne Twister 简介

• 迄今为止最好的随机数发生器之一
• 被标准C++、 Python、 MATLAB 等流行编程工具采用
• MT19937-32和MT19937-64的使用最为广泛
• 周期为2¹⁹⁹³⁷-1

不了解移位寄存器产生随机数的小伙伴可以去看这几个博客：
伪随机数生成算法-梅森旋转（Mersenne Twister/MT）算法介绍
伪随机数生成——梅森旋转（Mersenne Twister/MT）算法笔记

Mersenne Twister 的使用

1. Mersenne Twister 递推式
2. Twist
3. A的特殊结构使得xA的计算得到简化：
   1. xA = x >> 1 ; (x0 = 0)
   2. xA = (x >> 1) xor a ; (x0 = 1)
4. Temper: 为了改善均匀性引入的计算：
   1. y = x xor (x >> u)
   2. y = y xor ((y << s) & b)
   3. y = y xor ((y << t) & c)
   4. z = z xor (y >> l)

Mersenne Twister 相关代码

1.伪代码

//创建一个长度为624的数组来存储发生器的状态
int[0..623] MT int index = 0  
//用一个种子初始化发生器
//function initialize_generator(int seed) 
{i := 0  MT[0] := seed for i from 1 to 623 { // 遍历剩下的每个元素         MT[i] := last 32 bits of(1812433253 * (MT[i-1] xor (right shift by 30 bits(MT[i-1]))) + i) // 1812433253 == 0x6c078965     } 
}  // Extract a tempered pseudorandom number based on the index-th value,
// calling generate_numbers() every 624 numbers function extract_number() 
{     if index == 0 {    generate_numbers()    }      int y := MT[index]     y := y xor (right shift by 11 bits(y))    y := y xor (left shift by 7 bits(y) and (2636928640)) // 2636928640 == 0x9d2c5680    y := y xor (left shift by 15 bits(y) and (4022730752))// 4022730752 == 0xefc60000     y := y xor (right shift by 18 bits(y))     index := (index + 1) mod 624     return y 
}  // Generate an array of 624 untempered numbers
function generate_numbers() 
{for i from 0 to 623 {int y := (MT[i] & 0x80000000)                       // bit 31 (32nd bit) of MT[i]                     + (MT[(i+1) mod 624] & 0x7fffffff)   // bits 0-30 (first 31 bits) of MT[...]         MT[i] := MT[(i + 397) mod 624] xor (right shift by 1 bit(y))        if (y mod 2) != 0 { // y is odd             MT[i] := MT[i] xor (2567483615) // 2567483615 == 0x9908b0df         }     }}

2. Python 代码

#梅森旋转算法   -xlxw
#参考:mersenne twister from wikipedia#import
from time import time
import numpy as np#var
index = 624
MT = [0]*index
# MT[0] ->seeddef inter(t):return(0xFFFFFFFF & t) #取最后32位->tdef twister():global indexfor i in range(624):y = inter((MT[i] & 0x80000000) +(MT[(i + 1) % 624] & 0x7fffffff))MT[i] = MT[(i + 397) % 624] ^ y >> 1if y % 2 != 0:MT[i] = MT[i] ^ 0x9908b0dfindex = 0def exnum():global indexif index >= 624:twister()y = MT[index]y = y ^ y >> 11y = y ^ y << 7 & 2636928640y = y ^ y << 15 & 4022730752y = y ^ y >> 18index = index + 1return inter(y)def mainset(seed):MT[0] = seed    #seedfor i in range(1,624):MT[i] = inter(1812433253 * (MT[i - 1] ^ MT[i - 1] >> 30) + i)return exnum()def main():br = input("请输入随机数产生的范围(用,隔开):")mi = eval(br.split(',')[0])ma = eval(br.split(',')[1])    so = mainset(int(time())) / (2**32-1)rd = mi + int((ma-mi)*so)print("产生的随机整数为：",rd)
main()

3. Matlab 代码
根据Python就很容易可以写出MATLAB实现了

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！