From f459b195f4e979247e79fdbffe93ca5fee1f466d Mon Sep 17 00:00:00 2001
From: mhjensen <morten.hjorth-jensen@fys.uio.no>
Date: Thu, 17 May 2018 22:24:48 -0400
Subject: [PATCH] Added new codes

---
 doc/Programs/Finance/agents.py                |   7 +-
 doc/Programs/RandomWalks/OneDimParticle.py    |  45 ++++++++
 doc/Programs/RandomWalks/randomwalkfitting.py | 101 ++++++++++++++++++
 3 files changed, 151 insertions(+), 2 deletions(-)
 create mode 100644 doc/Programs/RandomWalks/OneDimParticle.py
 create mode 100644 doc/Programs/RandomWalks/randomwalkfitting.py

diff --git a/doc/Programs/Finance/agents.py b/doc/Programs/Finance/agents.py
index fad71e19b..8f5916a1d 100644
--- a/doc/Programs/Finance/agents.py
+++ b/doc/Programs/Finance/agents.py
@@ -1,3 +1,6 @@
+# Simulation of financial transations with or without saving/taxation on transaction
+# If lambda =0.0, no saving/taxation
+# See Patriarca et al http://www.sciencedirect.com/science/article/pii/S0378437104004327
 #!/usr/bin/env python
 import numpy as np
 import matplotlib.mlab as mlab
@@ -11,7 +14,7 @@ Agents  = 500
 MCcounts = 1000
 Transactions = 100000
 startMoney = 1.0
-Lambda = 0.0
+Lambda = 0.2
 FinancialAgents = startMoney*np.ones(Agents)
 for i in range (1, MCcounts, 1):
     for j in range (1, Transactions, 1):
@@ -30,6 +33,6 @@ n, bins, patches = plt.hist(FinancialAgents, 50, facecolor='green')
 plt.xlabel('$x$')
 plt.ylabel('Distribution of wealth')
 plt.title(r'Money')
-plt.axis([0, 10, 0, 500])
+plt.axis([0, 10, 0, 100])
 plt.grid(True)
 plt.show()
diff --git a/doc/Programs/RandomWalks/OneDimParticle.py b/doc/Programs/RandomWalks/OneDimParticle.py
new file mode 100644
index 000000000..ec77cdb42
--- /dev/null
+++ b/doc/Programs/RandomWalks/OneDimParticle.py
@@ -0,0 +1,45 @@
+# Program to test the Metropolis algorithm with one particle at given temp in
+# one dimension
+#!/usr/bin/env python
+import numpy as np
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+import random
+from math import sqrt, exp, log
+# initialize the rng with a seed
+random.seed()
+# Hard coding of input parameters
+MCcycles = 100000
+Temperature = 2.0
+beta = 1./Temperature
+InitialVelocity = -2.0
+CurrentVelocity = InitialVelocity
+Energy = 0.5*InitialVelocity*InitialVelocity
+VelocityRange = 10*sqrt(Temperature)
+VelocityStep = 2*VelocityRange/10.
+AverageEnergy = Energy
+AverageEnergy2 = Energy*Energy
+VelocityValues = np.zeros(MCcycles)
+# The Monte Carlo sampling with Metropolis starts here
+for i in range (1, MCcycles, 1):
+    TrialVelocity = CurrentVelocity + (2.0*random.random() - 1.0)*VelocityStep
+    EnergyChange = 0.5*(TrialVelocity*TrialVelocity -CurrentVelocity*CurrentVelocity);
+    if random.random() <= exp(-beta*EnergyChange):
+        CurrentVelocity = TrialVelocity
+        Energy += EnergyChange
+        VelocityValues[i] = CurrentVelocity
+    AverageEnergy += Energy
+    AverageEnergy2 += Energy*Energy
+#Final averages
+AverageEnergy = AverageEnergy/MCcycles
+AverageEnergy2 = AverageEnergy2/MCcycles
+Variance = AverageEnergy2 - AverageEnergy*AverageEnergy
+print(AverageEnergy, Variance)
+n, bins, patches = plt.hist(VelocityValues, 400, facecolor='green')
+
+plt.xlabel('$v$')
+plt.ylabel('Velocity distribution P(v)')
+plt.title(r'Velocity histogram at $k_BT=2$')
+plt.axis([-5, 5, 0, 600])
+plt.grid(True)
+plt.show()
diff --git a/doc/Programs/RandomWalks/randomwalkfitting.py b/doc/Programs/RandomWalks/randomwalkfitting.py
new file mode 100644
index 000000000..bd4a62239
--- /dev/null
+++ b/doc/Programs/RandomWalks/randomwalkfitting.py
@@ -0,0 +1,101 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.linear_model import LinearRegression
+
+steps=250
+
+distance=0
+x=0
+distance_list=[]
+steps_list=[]
+while x<steps:
+    distance+=np.random.randint(-1,2)
+    distance_list.append(distance)
+    x+=1
+    steps_list.append(x)
+plt.plot(steps_list,distance_list, color='green', label="Random Walk Data")
+
+steps_list=np.asarray(steps_list)
+distance_list=np.asarray(distance_list)
+
+X=steps_list[:,np.newaxis]
+
+#Polynomial fits
+
+#Degree 2
+poly_features=PolynomialFeatures(degree=2, include_bias=False)
+X_poly=poly_features.fit_transform(X)
+
+lin_reg=LinearRegression()
+poly_fit=lin_reg.fit(X_poly,distance_list)
+b=lin_reg.coef_
+c=lin_reg.intercept_
+print ("2nd degree coefficients:")
+print ("zero power: ",c)
+print ("first power: ", b[0])
+print ("second power: ",b[1])
+
+z = np.arange(0, steps, .01)
+z_mod=b[1]*z**2+b[0]*z+c
+
+fit_mod=b[1]*X**2+b[0]*X+c
+plt.plot(z, z_mod, color='r', label="2nd Degree Fit")
+plt.title("Polynomial Regression")
+
+plt.xlabel("Steps")
+plt.ylabel("Distance")
+
+#Degree 10
+poly_features10=PolynomialFeatures(degree=10, include_bias=False)
+X_poly10=poly_features10.fit_transform(X)
+
+poly_fit10=lin_reg.fit(X_poly10,distance_list)
+
+y_plot=poly_fit10.predict(X_poly10)
+plt.plot(X, y_plot, color='black', label="10th Degree Fit")
+
+plt.legend()
+plt.show()
+
+
+#Decision Tree Regression
+from sklearn.tree import DecisionTreeRegressor
+regr_1=DecisionTreeRegressor(max_depth=2)
+regr_2=DecisionTreeRegressor(max_depth=5)
+regr_3=DecisionTreeRegressor(max_depth=7)
+regr_1.fit(X, distance_list)
+regr_2.fit(X, distance_list)
+regr_3.fit(X, distance_list)
+
+X_test = np.arange(0.0, steps, 0.01)[:, np.newaxis]
+y_1 = regr_1.predict(X_test)
+y_2 = regr_2.predict(X_test)
+y_3=regr_3.predict(X_test)
+
+# Plot the results
+plt.figure()
+plt.scatter(X, distance_list, s=2.5, c="black", label="data")
+plt.plot(X_test, y_1, color="red",
+         label="max_depth=2", linewidth=2)
+plt.plot(X_test, y_2, color="green", label="max_depth=5", linewidth=2)
+plt.plot(X_test, y_3, color="m", label="max_depth=7", linewidth=2)
+
+plt.xlabel("Data")
+plt.ylabel("Darget")
+plt.title("Decision Tree Regression")
+plt.legend()
+plt.show()
+
+
+"""new_dist=distance_list[-1]
+step_max=2500
+new_x=steps
+new_dist_list=[]
+new_steps_list=np.arange(steps,step_max)
+while new_x>=steps and new_x<step_max:
+    dist_prediction=clf.predict(new_x)
+    new_dist_list.append(dist_prediction)
+    new_x+=1
+plt.plot(new_steps_list,new_dist_list, color='red')
+plt.show()"""