Skip to content

Commit

Permalink
Newly installed git was triggered by newlines
Browse files Browse the repository at this point in the history
  • Loading branch information
weka511 committed Jun 2, 2022
1 parent c607aa0 commit 1895b66
Show file tree
Hide file tree
Showing 251 changed files with 23,870 additions and 23,860 deletions.
66 changes: 33 additions & 33 deletions book/binomialconvolution.py
Original file line number Diff line number Diff line change
@@ -1,34 +1,34 @@
# binomialconvolution.py

# Copyright (C) 2015,2018 Greenweaves Software Pty Ltd

# This is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This software is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>

import math,matplotlib.pyplot as plt

# Algorithm 1.25 from Krauth
def binomial_convolution(theta,pi):
def calculate_row(pi):
return [theta*x + (1-theta)*y for (x,y) in zip(pi[:-1],pi[1:])]
return calculate_row([0]+pi+[0])


if __name__=="__main__":
pi=[1]
theta=math.pi/4
for i in range(8):
pi=binomial_convolution(theta,pi)
plt.plot(pi)
plt.show()
# binomialconvolution.py

# Copyright (C) 2015,2018 Greenweaves Software Pty Ltd

# This is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This software is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>

import math,matplotlib.pyplot as plt

# Algorithm 1.25 from Krauth
def binomial_convolution(theta,pi):
def calculate_row(pi):
return [theta*x + (1-theta)*y for (x,y) in zip(pi[:-1],pi[1:])]
return calculate_row([0]+pi+[0])


if __name__=="__main__":
pi=[1]
theta=math.pi/4
for i in range(8):
pi=binomial_convolution(theta,pi)
plt.plot(pi)
plt.show()

62 changes: 31 additions & 31 deletions book/boxmuller.py
Original file line number Diff line number Diff line change
@@ -1,31 +1,31 @@
# Copyright (C) 2015-2012 Greenweaves Software Limited

# This is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This software is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.

''' Test for Box Muller algorithm'''

from matplotlib.pyplot import hist, show, title, xlabel, ylabel
from smac import BoxMuller

if __name__=="__main__":

gauss = BoxMuller()

hist([gauss.gauss() for _ in range(1000000)],
bins = 200,
density = True)
title("Gaussian Histogram")
xlabel("Value")
ylabel("Frequency")
show()
# Copyright (C) 2015-2012 Greenweaves Software Limited

# This is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This software is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.

''' Test for Box Muller algorithm'''

from matplotlib.pyplot import hist, show, title, xlabel, ylabel
from smac import BoxMuller

if __name__=="__main__":

gauss = BoxMuller()

hist([gauss.gauss() for _ in range(1000000)],
bins = 200,
density = True)
title("Gaussian Histogram")
xlabel("Value")
ylabel("Frequency")
show()
70 changes: 35 additions & 35 deletions book/direct-gamma-zeta.py
Original file line number Diff line number Diff line change
@@ -1,36 +1,36 @@
# Copyright (C) 2015 Greenweaves Software Pty Ltd

# This is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>

import random,math

def direct_gamma_zeta(gamma,zeta,n):
sigma=0
sigma2=0
for i in range(n):
x=0
while x==0:
x=random.random()
x1=x**(1/(1+zeta))
x2=x1**(gamma-zeta)
sigma+=x2
sigma2+=(x2*x2)
mean=sigma/n
return (mean,math.sqrt(sigma2/n-mean*mean)/math.sqrt(n))

def run(zeta,n):
print(zeta)
for gamma in [2.0,1.0,0.0,-0.1,-0.4,-0.8]:
(s,t)=direct_gamma_zeta(gamma,zeta,n)
print (gamma, s-t,s+t, (zeta+1)/(gamma+1))

if __name__=="__main__":
for zeta in [0.0, -0.1, -0.7]:
# Copyright (C) 2015 Greenweaves Software Pty Ltd

# This is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>

import random,math

def direct_gamma_zeta(gamma,zeta,n):
sigma=0
sigma2=0
for i in range(n):
x=0
while x==0:
x=random.random()
x1=x**(1/(1+zeta))
x2=x1**(gamma-zeta)
sigma+=x2
sigma2+=(x2*x2)
mean=sigma/n
return (mean,math.sqrt(sigma2/n-mean*mean)/math.sqrt(n))

def run(zeta,n):
print(zeta)
for gamma in [2.0,1.0,0.0,-0.1,-0.4,-0.8]:
(s,t)=direct_gamma_zeta(gamma,zeta,n)
print (gamma, s-t,s+t, (zeta+1)/(gamma+1))

if __name__=="__main__":
for zeta in [0.0, -0.1, -0.7]:
run(zeta,1000)
126 changes: 63 additions & 63 deletions book/direct-gamma.py
Original file line number Diff line number Diff line change
@@ -1,63 +1,63 @@
# MIT License
#
# Copyright (c) 2018 Simon Crase
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Implement Algorithm 1.29, subtract mean value for each sample, and generate
# histograms of the average of N samples and the rescaled averages.

import random,pylab,math,numpy as np
from scipy.stats import gaussian_kde

def direct_gamma(gamma,generator=False,K=10000,N=1000000):
sigma = 0
for i in range(N):
sigma += random.random()**gamma
return sigma/N


if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser(description='Evaluate integral from Krauth Section 1.4.')
parser.add_argument('steps', metavar='M', type=int, nargs=1,help='Number of steps for integral')
parser.add_argument('--N', metavar='N', type=int, nargs='+',default=[1,10,100,1000,10000],help='Number of steps for integral')
parser.add_argument('--gamma',metavar='gamma',type=float,nargs=1,default=-0.8,help='exponent')
args = parser.parse_args()
M = args.steps[0]
gamma = args.gamma
bins=[0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,100000000000]

for N in args.N:
def scaled(n):
return (direct_gamma(gamma,N=n)-5)/(N**(-1-gamma))
data = [direct_gamma(gamma,N=N) for m in range(M)]
density = gaussian_kde(data)
y,binEdges=np.histogram(data,normed=True,bins=bins)
bincenters = 0.5*(binEdges[1:]+binEdges[:-1])
pylab.plot(bincenters,density(bincenters),label='N={0}'.format(N))
pylab.xlim(1,10)
pylab.xlabel(r'$\Sigma/N$')
pylab.ylabel(r'$\pi(\Sigma/N$)')
pylab.title(r'$Average\ {0}\ iterations$'.format(M))
pylab.legend()
pylab.savefig('histogram.png')
pylab.show()


# MIT License
#
# Copyright (c) 2018 Simon Crase
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Implement Algorithm 1.29, subtract mean value for each sample, and generate
# histograms of the average of N samples and the rescaled averages.

import random,pylab,math,numpy as np
from scipy.stats import gaussian_kde

def direct_gamma(gamma,generator=False,K=10000,N=1000000):
sigma = 0
for i in range(N):
sigma += random.random()**gamma
return sigma/N


if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser(description='Evaluate integral from Krauth Section 1.4.')
parser.add_argument('steps', metavar='M', type=int, nargs=1,help='Number of steps for integral')
parser.add_argument('--N', metavar='N', type=int, nargs='+',default=[1,10,100,1000,10000],help='Number of steps for integral')
parser.add_argument('--gamma',metavar='gamma',type=float,nargs=1,default=-0.8,help='exponent')
args = parser.parse_args()
M = args.steps[0]
gamma = args.gamma
bins=[0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,100000000000]

for N in args.N:
def scaled(n):
return (direct_gamma(gamma,N=n)-5)/(N**(-1-gamma))
data = [direct_gamma(gamma,N=N) for m in range(M)]
density = gaussian_kde(data)
y,binEdges=np.histogram(data,normed=True,bins=bins)
bincenters = 0.5*(binEdges[1:]+binEdges[:-1])
pylab.plot(bincenters,density(bincenters),label='N={0}'.format(N))
pylab.xlim(1,10)
pylab.xlabel(r'$\Sigma/N$')
pylab.ylabel(r'$\pi(\Sigma/N$)')
pylab.title(r'$Average\ {0}\ iterations$'.format(M))
pylab.legend()
pylab.savefig('histogram.png')
pylab.show()

Loading

0 comments on commit 1895b66

Please sign in to comment.