擬似乱数の生成 by C

SetDirectory["/Users/takasu/home/化学生物環境学科の仕事/講義/2022 (R4) 年度/環境科学計算機実験 2022R4/疑似乱数生成/simulation/"]

一様乱数 Uniformly distributed 1

data=ReadList["data-1d",Real];
Length[data]

ListPlot[data]

Take[data,10]

Min[data]
Max[data]

gHisto=Histogram[data,{0,1,0.05},"PDF"]

myFunc=PDF[UniformDistribution[{0,1}],x]

gPDF=Plot[myFunc,{x,-1,2},PlotRange->All,PlotStyle->Thick]

Show[gHisto,gPDF]

一様乱数 Uniformly distributed 2

data=ReadList["data-2d",{Real,Real}];
Length[data]

Take[data,10]

ListPlot[data,AspectRatio->1]

一様乱数 Uniformly distributed 3

data=ReadList["data-3d",{Real,Real,Real}];
Length[data]

Take[data,10]

seqPoint=Map[Point,data];
Show[Graphics3D[seqPoint]]

ListPointPlot3D[data,BoxRatios->{1,1,1}]

指数乱数 Exponentially distributed

data=ReadList["data-exp",Real];
Length[data]

Take[data,10]

Min[data]
Max[data]

gHisto=Histogram[data,{0,5,0.1},"PDF"]

myFunc=PDF[ExponentialDistribution[2],x]

gPDF=Plot[myFunc,{x,0,5},PlotRange->All]

Show[gHisto,gPDF]

正規乱数 Normally distributed

data=ReadList["data-gas",Real];
Length[data]

Take[data,10]

Min[data]
Max[data]

gHisto=Histogram[data,{-4,4,0.2},"PDF"]

myFunc=PDF[NormalDistribution[0,1],x]

gPDF=Plot[myFunc,{x,-4,4}]

Show[gHisto,gPDF]

疑似乱数の生成 by Mathematica

一様乱数（0~ 1 の実数）

Random[]

(* U[0,1) を 10000個生成する *)
len = 10000;

seq={};
Do[
AppendTo[seq,Random[]],{i,1,len}
]

(* はじめの100個を表示 *)
Take[seq,100]

Histogram[seq,Automatic]

(* ヒストグラムを確率密度関数として描く *)
gHist=Histogram[seq,Automatic,"PDF"]

(* 一様乱数の確率密度関数を描く *)
gPdf=Plot[Evaluate[PDF[UniformDistribution[{0,1}],k]],{k,-0.2,1.2},PlotRange->All]

Show[gHist,gPdf]

確率変数の変換 0

Mathematicaを用いた0 - 1 の一様乱数の生成

Random[]

(* U[0,1) を 10000個生成する *)
len = 10000;

seq={};
Do[
AppendTo[seq,Random[]],{i,1,len}
]

Take[seq,10]

(* seqの各要素を変換する *)
seq2=Sqrt[seq]

Take[seq2,10]

Max[seq2]

Min[seq2]

gHist=Histogram[seq2,Automatic,"PDF"]

確率変数の変換 1 Exponential distribution

Mathematicaを用いた0 - 1 の一様乱数の生成

Random[]

(* U[0,1) を 10000個生成する *)
len = 10000;

seq={};
Do[
AppendTo[seq,Random[]],{i,1,len}
]

Take[seq,10]

(* 対数をとって符号を変える *)
seq2=-Log[seq]/2

Take[seq2,10]

Max[seq2]

Min[seq2]

gHist=Histogram[seq2,Automatic,"PDF"]

gPdf=Plot[Evaluate[PDF[ExponentialDistribution[2],k]],{k,0,10},PlotRange->All]

Show[gHist,gPdf]

確率変数の変換 2 Normal distribution, Box Muller 法

(* U[0,1) を 10000個生成する *)
len = 10000;

seq1={};
Do[
AppendTo[seq1,Random[]],{i,1,len}
]
seq2={};
Do[
AppendTo[seq2,Random[]],{i,1,len}
]

(* 平方根をとる *)
seq3=Sqrt[-2Log[seq1]]Cos[2Pi seq2];

Take[seq3,10]

Max[seq3]

Min[seq3]

gHist=Histogram[seq3,{-4,4,0.2},"PDF"]

x=.

gPdf=Plot[Evaluate[PDF[NormalDistribution[0,1],x]],{x,-3,3},PlotRange->All]

Show[gHist,gPdf]

pdfGauss=Exp[-(x-m)^2/2/sd^2]/Sqrt[2Pi]/sd/.{m->0,sd->1}

Plot[pdfGauss,{x,-3,3}]

Integrate[pdfGauss,{x,-Infinity,Infinity}]

seq4=seq3 2+1

g=Histogram[seq4,Automatic,"PDF"]

gPdf=Plot[Evaluate[PDF[NormalDistribution[1,2],x]],{x,-10,10},PlotRange->All]

Show[g,gPdf]

確率変数の変換 3

(* U[0,1) を 10000個生成する *)
len = 10000;

seq={};
Do[
AppendTo[seq,Random[]],{i,1,len}
]

Take[seq,10]

(* 平方根をとる *)
seq2=Sqrt[seq];

Take[seq2,10]

gHist=Histogram[seq2,Automatic,"PDF"]

gPdf=Plot[2k,{k,0,1},PlotRange->All]

Show[gHist,gPdf]

ベルヌイ試行に従う乱数（離散確率変数）の生成

確率 P (0 <= P <= 1) で起こる事象 A を考える（ベルヌイ試行）

bernoulli[p_]:=Module[{rv, result},
rv=Random[];
If[rv ≤p, result=1,result=0];
Return[result];
]

bernoulli[0.5]

(* パラメータ p のベルヌイ試行に従う乱数を 10000 個生成 *)
len = 10000;
p=1/2;

seq={};
Do[
AppendTo[seq,bernoulli[p]],{i,1,len}
]

(* はじめの100個を表示 *)
Take[seq,10]

Count[seq,0]

Count[seq,1]

(* ヒストグラムを確率密度関数として描く *)
gHist=Histogram[seq,Automatic,"PDF"]

二項分布に従う乱数（離散確率変数）の生成

Random[BinomialDistribution[10,0.5]]

(* パラメータ n, p の二項分布に従う乱数を 10000 個生成 MathematicaのBinomialDistributionを用いる *)
len = 10000;
n=20;
p=1/2;

seq={};
Do[
AppendTo[seq,Random[BinomialDistribution[n,p]]],{i,1,len}
]

(* はじめの100個を表示 *)
Take[seq,10]

(* ヒストグラムを確率密度関数として描く *)
gHist=Histogram[seq,n,"PDF"]

Table[Evaluate[PDF[BinomialDistribution[n,p],k]],{k,0,n}]

(* 確率密度関数を描く *)
gPdf=DiscretePlot[Evaluate[PDF[BinomialDistribution[n,p],k]],{k,n},PlotRange->All,PlotMarkers->Automatic]

Show[gHist,gPdf]

(* 一様乱数 Random[] を用いて定義する二項分布乱数 *)

binomial[n_,p_]:=Module[{i,count=0},
Do[
If[Random[]<p,count++],{i,1,n}
];
Return[count];
]

binomial[10,0.5]

seq2=Table[binomial[10,0.5],10000];

n

(* ヒストグラムを確率密度関数として描く *)
gHist2=Histogram[seq2,10,"PDF"]

(* 確率密度関数を描く *)
p = 0.5;
n=10;
gPdf2=DiscretePlot[Evaluate[PDF[BinomialDistribution[n,p],k]],{k,n},PlotRange->All,PlotMarkers->Automatic]

Show[gHist2,gPdf2]

ポアソン Poisson 分布に従う乱数（離散確率変数）の生成 　20210421追加

Random[PoissonDistribution[3]]

(* パラメータ lambda のPoisson分布に従う乱数を 10000 個生成 MathematicaのPoissonDistributionを用いる *)
len = 10000;

lambda=3;

seq={};
Do[
AppendTo[seq,Random[PoissonDistribution[lambda]]],{i,1,len}
]

(* はじめの100個を表示 *)
Take[seq,10]

(* ヒストグラムを確率密度関数として描く *)
gHist=Histogram[seq,{0,10,1},"PDF"]

平均 lambda のポアソン分布の理論値

Table[Evaluate[PDF[PoissonDistribution[lambda],k]],{k,0,10}]

(* 確率密度関数を描く *)
gPdf=DiscretePlot[Evaluate[PDF[PoissonDistribution[lambda],k]],{k,10},PlotRange->All,PlotMarkers->Automatic]

Show[gHist,gPdf]

モンテカルロ積分

円周率の近似値として半径1の円の面積を求める

2 Integrate[Sqrt[1-x^2],{x,-1,1}]

N[%]

countTrial=0;
countIn=0;
While[countTrial<1000000,
x=Random[Real,{-1,1}];
y=Random[Real,{-1,1}];
If[x^2+y^2<1,countIn++];
countTrial++;
];

countIn

pi=N[4 countIn/countTrial]

モンテカルロシミュレーション

Random walk on lattice by Mathematica

numIndivs=100;
indivs=Table[{0,0},{i,1,numIndivs}];

tStep=100;
data={};

Do[
AppendTo[data,indivs];

Do[
dx=0;
dy=0;

r=RandomReal[{0,1}];
If[r<1/4,
dx=1,
If[r<1/2,
dx=-1,
If[r<3/4,
dy=1,dy=-1]
]
];

indivs[[i]]+={dx,dy},{i,1,numIndivs}

]

,{t,0,tStep}
]

Take[Last[data],10]

glist={};
Do[
g=
ListPlot[data[[i]],PlotRange->{{-30,30},{-30,30}},AspectRatio->1,PlotStyle->{PointSize[0.01]}];
AppendTo[glist,g]
,{i,1,Length[data]}
]

ListAnimate[glist]

Export["RandomWalkLattice.mov",glist]

?Export

maxDistance=Map[Max,data]

ListPlot[maxDistance,Joined->True,AxesLabel->{"t","Max distance"},PlotStyle->Thickness[0.005]]

Random walk by Mathematica

numIndivs=100;
indivs=Table[{0,0},{i,1,numIndivs}];

tStep=100;

d=1;
data={};
Do[
AppendTo[data,indivs];

Do[
angle=RandomReal[{0,2Pi}];
{dx,dy}={d Cos[angle],d Sin[angle]};
indivs[[i]]+={dx,dy},{i,1,numIndivs}
]
,{t,0,tStep}
]

Take[Last[data],10]

glist={};
Do[
g=
ListPlot[data[[i]],PlotRange->{{-30,30},{-30,30}},AspectRatio->1,PlotStyle->PointSize[0.01]];
AppendTo[glist,g],{i,1,Length[data]}
]

ListAnimate[glist]

maxDistance=Map[Max,data]

ListPlot[maxDistance,Joined->True,AxesLabel->{"t","Max distance"},PlotStyle->Thickness[0.005]]

Random walk by C

SetDirectory["~/"]

data=ReadList["data_rw",Real,RecordLists->True];
data2=Map[Partition[#,2]&,data];
Length[data2]
{Max[data2],Min[data2]}

glist={};
Do[
g=ListPlot[data2[[i]],PlotRange->{{-30,30},{-30,30}},AspectRatio->1,PlotStyle->{PointSize[0.01]}];
AppendTo[glist,g],{i,1,Length[data2]}
]

ListAnimate[glist]

glist

glist[[10]]