
# import packages that provides certain function
# install a package: pip install packageName
# uninstall a package: pip uninstall packageName

import math 
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.integrate import quad
from sympy import symbols, diff, sin, cos


# The text behind the symbol # is a comment in Python
1 + 4

5


# math.pi is circumference ratio which is approximately 3.1415926 
2 + math.pi/4-0.8

1.9853981633974482


x = 1 # "=" is the assignment symbol: it assigns 1 to x.
y = 2
z = 4

# This statement is equivalent to 2*1^2-4 = -2
t = 2*x**y-z
t

-2


u = 2
v = 3
u + v

5


math.sin(u*v) # u*v = 6 is considered radian

-0.27941549819892586


# This statement creates function where it is a block of code
# that runs when it is called.
# It may pass in data which are called parameters

def square(x):
    return x*x
square(4)

16


def fctn(x, y, z):
    return x + y - z/2

fctn(1, 2, 3)

1.5


# Plot the curve of y = sin(x) from -pi to 2*pi
x = np.linspace(-np.pi, 2*np.pi, 1000)
plt.plot(x, np.sin(x), color='k')
plt.xlabel('x')
plt.ylabel('$\sin(x)$')
plt.show()
plt.clf()

<Figure size 640x480 with 0 Axes>


# Define a function
def square(x):
    return x*x

# Plot the defined function
plt.plot(np.linspace(-3, 2, 100), square(np.linspace(-3, 2, 100)), color='k')
plt.xlabel('x')
plt.ylabel('$x^{2}$')
plt.show()
plt.clf()

<Figure size 640x480 with 0 Axes>


from mpl_toolkits import mplot3d
# We import the mplot3d from mpl_toolkits, another python package necessary for plotting in 3D.

def f(x,y):
    return 1 - x**2 - y**2


# Plot a 3D surface with perspective angle 310 deg
x = y = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x, y)
Z = f(X,Y)
fig = plt.figure(figsize=(16,12))
ax = plt.axes(projection='3d')
ax.contour3D(X, Y, Z, 200, cmap='binary')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
plt.show()


x = symbols('x') #creates symbolic 'x' variable

diff(x**2, x) # Take derivative of x^2 which is 2x


fx = x**2 # define a symbolic function with variable x

diff(fx, x) # differentiate the function with respect to x


# Change the expression and use the same derivative command
fx = x**2 * sin(x)
diff(fx, x)


#create 3 symbolic variables to create a multi-variable function
x, y, z = symbols('x y z')
fxyz = x**2 + y**2 + z**2

fxyz # This gives the expression of the function in x, y and z


diff(fxyz, x) # gives the partial derivative with respect to x


diff(fxyz, y) # gives the partial derivative with respect to y


def square(x):
    return x**2

# Numerically integrating x^2 from 0 to 1
result, abs_error = quad(square, 0, 1)
print("Integration result:", result)
print("Absolute error:", abs_error)

Integration result: 0.33333333333333337
Absolute error: 3.700743415417189e-15


# Integrates cos(x) from 0 to pi/2
result, abs_error = quad(cos, 0, math.pi/2)
print("Integration result:", result)
print("Absolute error:", abs_error)

Integration result: 0.9999999999999999
Absolute error: 1.1102230246251564e-14


# Enter data inside np.array() to create a 5x1 column vector
np.array([1, 6, 3, math.pi, -3])

array([ 1.        ,  6.        ,  3.        ,  3.14159265, -3.        ])


# Generate a sequence from 2 to 6
np.arange(2, 7) #not upper boundary inclusive

array([2, 3, 4, 5, 6])


# Generate a sequence from 1 to 10 with an increment of 2
np.arange(1, 11, 2)

array([1, 3, 5, 7, 9])


x = np.array([[1, -1, 1, -1]])
# 1 is added to each element of x
x + 1

array([[2, 0, 2, 0]])


# 2 multiplies each element of x
2 * x

array([[ 2, -2,  2, -2]])


# Each element of x is divided by 2
x / 2

array([[ 0.5, -0.5,  0.5, -0.5]])


y = np.array([np.arange(1, 5)])
x * y # * multiplies each pair of elements

array([[ 1, -2,  3, -4]])

x

array([[ 1, -1,  1, -1]])


# Transpose of x matrix
x.T

array([[ 1],
       [-1],
       [ 1],
       [-1]])


# 4x1 matrix times a 1x4 matrix yields a 4x4 matrix
np.dot(x,y.T)

array([[-2]])


# Matrix multiplication: 1x4 matrix times a 4x1 matrix
# equivalent to a dot product
np.dot(x.T, y)

array([[ 1,  2,  3,  4],
       [-1, -2, -3, -4],
       [ 1,  2,  3,  4],
       [-1, -2, -3, -4]])


mx = x.reshape(2,-1)
mx

array([[ 1, -1],
       [ 1, -1]])


my = y.reshape(2,-1)
my

array([[1, 2],
       [3, 4]])


# Dimensions of a matrix
np.shape(my)

(2, 2)


# Convert a matrix to a vector
my.flatten()

array([1, 2, 3, 4])


# Multiplication between each pair of elements
mx * my

array([[ 1, -2],
       [ 3, -4]])


# Division between each pair of elements
mx / my

array([[ 1.        , -0.5       ],
       [ 0.33333333, -0.25      ]])


mx - 2 * my

array([[-1, -5],
       [-5, -9]])


# Multiplication of two matrices
np.dot(my, mx)

array([[ 3, -3],
       [ 7, -7]])


# Determinant of a square matrix
np.linalg.det(my)

-2.0000000000000004


# Find inverse of a matrix
myinv = np.linalg.inv(my)
myinv

array([[-2. ,  1. ],
       [ 1.5, -0.5]])


# Verifies the inverse of a matrix and results in identity matrix
np.dot(myinv, my)

array([[1.00000000e+00, 0.00000000e+00],
       [1.11022302e-16, 1.00000000e+00]])


# outputs the diagonal vector of a matrix
np.diag(my)

array([1, 4])


# np.linalg.eig computes the eigenvalues and eigenvectors of a matrix.
# Each are returned as separate arrays, so we store each in a variable to display them individually
eigenValues, eigenVectors = np.linalg.eig(my)
print(f"\nEigenvalues:\n{eigenValues}")
print(f"\nEigenvectors:\n{eigenVectors}")

Eigenvalues:
[-0.37228132  5.37228132]

Eigenvectors:
[[-0.82456484 -0.41597356]
 [ 0.56576746 -0.90937671]]


# np.linalg.svd computes the singular value decomposition of a matrix.
# This one returns multiple arrays, so we store each in it's own variable and display them in turn.
U, S, V = np.linalg.svd(my)
print(f"\nLeft-singular matrix of my:\n{U}")
print(f"\nSingular values of my:\n{S}")
print(f"\nRight-singular matrix of my:\n{V}")

Left-singular matrix of my:
[[-0.40455358 -0.9145143 ]
 [-0.9145143   0.40455358]]

Singular values of my:
[5.4649857  0.36596619]

Right-singular matrix of my:
[[-0.57604844 -0.81741556]
 [ 0.81741556 -0.57604844]]


# generates 8 random numbers in (-infty, infty)
# with normal distribution with default mean = 0, and default sd = 1
x = np.random.randn(8)
x

array([0.59046501, 1.34782775, 0.62694335, 0.70552022, 1.89571115,
       1.11598687, 0.18258704, 0.03792304])


# mean
np.mean(x)


# variance
np.var(x)


# standard deviation
np.std(x)


# median
np.median(x)


# range
np.min(x), np.max(x)


plt.boxplot(x, patch_artist=True,
            boxprops = dict(facecolor = 'lightgrey'),
            medianprops = dict(color = "black", linewidth = 1.5),
            whiskerprops = dict(linestyle=':',linewidth=1.0,
                                color='black'))
plt.xticks([])
plt.show()


w = np.random.normal(size=1000)
z = np.random.normal(loc=10, scale=5, size=1000)


pd.DataFrame(np.random.randn(12)).describe()


# plot the histogram of 1000 random numbers
plt.hist(w, bins = 9, edgecolor = "black",color="lightgrey")
plt.xlabel('w')
plt.ylabel('Frequency')
plt.title('Histogram of w',fontweight='bold')
plt.show()


# Display working directory
import os
print(os.getcwd())

/Users/momtaza/Desktop/RMathModel


# Change working directory
os.chdir('/Users/momtaza/Desktop/RMathModel/data')


# Commands to input data files:

mydata = pd.read_csv("/Users/momtaza/Desktop/RMathModel/data/mydata.csv")
# read csv file named "mydata.csv"

mydata = pd.read_csv("/Users/momtaza/Desktop/RMathModel/data/mydata.txt", sep='\s+')
# read text file named "mydata.txt"

mydata = pd.read_excel("/Users/momtaza/Desktop/RMathModel/data/mydata.xlsx")
# read an excel file


import netCDF4 as nc

# Open a NetCDF file
ncin = nc.Dataset("air.mon.mean.nc", "r")  # 'r' is for read mode

# Read data "lon" from a NetCDF file
lon = ncin.variables['lon'][:]
print(lon)

# Close the NetCDF file after use
ncin.close()

[  0.    2.5   5.    7.5  10.   12.5  15.   17.5  20.   22.5  25.   27.5
  30.   32.5  35.   37.5  40.   42.5  45.   47.5  50.   52.5  55.   57.5
  60.   62.5  65.   67.5  70.   72.5  75.   77.5  80.   82.5  85.   87.5
  90.   92.5  95.   97.5 100.  102.5 105.  107.5 110.  112.5 115.  117.5
 120.  122.5 125.  127.5 130.  132.5 135.  137.5 140.  142.5 145.  147.5
 150.  152.5 155.  157.5 160.  162.5 165.  167.5 170.  172.5 175.  177.5
 180.  182.5 185.  187.5 190.  192.5 195.  197.5 200.  202.5 205.  207.5
 210.  212.5 215.  217.5 220.  222.5 225.  227.5 230.  232.5 235.  237.5
 240.  242.5 245.  247.5 250.  252.5 255.  257.5 260.  262.5 265.  267.5
 270.  272.5 275.  277.5 280.  282.5 285.  287.5 290.  292.5 295.  297.5
 300.  302.5 305.  307.5 310.  312.5 315.  317.5 320.  322.5 325.  327.5
 330.  332.5 335.  337.5 340.  342.5 345.  347.5 350.  352.5 355.  357.5]


import json

# Load/read JSON data from a file
with open("argo4903285_59.json", "r") as file:
    jd = json.load(file)

	0
count	12.000000
mean	-0.566646
std	1.035952
min	-3.105326
25%	-0.770634
50%	-0.593415
75%	0.059437
max	0.866987

Introduction to Modern Mathematical Modeling with R:

A User's Manual to Train Mathematical Consultants

A Cambridge University Press book

by SSP Shen

Chapter 2: Basics of R Programming¶

R as a Smart Calculator¶

Write a function in R¶

Plot with R¶

Symbolic Calculation for Calculus by R¶

Vectors and Matrices¶

Statistics¶

Python Tutorials¶