# coding: utf-8
# In[9]:
get_ipython().magic('matplotlib inline')
import matplotlib
import numpy as np
import matplotlib.pyplot as pt
# In[10]:
k = 2.4 #2pi/k is the frequency of sine wave
n = 100
x = np.linspace(0,10,n)
wave = np.sin(k*x) #create pure sine function or wave
noise = np.random.normal(scale = 1.0/10, size = len(x)) #normally distributed random numbers
y = wave + noise #add noise to curve
# In[11]:
#plot the curve
pt.plot(x,y,'o')
pt.xlim([-.5,10.5])
#pt.plot(x,np.sin(x))
#pt.plot(x,np.sin(2*x))
#pt.plot(x,np.sin(3*x))
pt.show()
# What are some properties of the function that give clues to what function it may be?
# In[23]:
# define residual of data and sine function as function of parameter k (sum of squares)
def r(k):
return sum((y - np.sin(k*x))**2)
# define derivative of residual with respect to k. (k is the unknown variable)
def rprime(k):
return -2*sum(x*np.cos(k*x)*(y-np.sin(k*x)))
# define 2nd derivative with respect to k
def r2prime(k):
# this gets really nasty so I'm going to define it in steps
s1 = (y - np.sin(k*x))*np.sin(k*x) + (np.cos(k*x))**2
s2 = s1*x**2
return 2*sum(s2)
# In[24]:
#Let's look at a plot of the residual
k = np.linspace(0,3,n)#it looks like the true k falls in this range
res = np.zeros(len(k))
for i in range(len(k)):
res[i] = r(k[i])
pt.plot(k,res)
pt.xlabel("value of parameter k")
pt.ylabel("residual")
pt.show()
# In[27]:
# there are a lot of local minima, so we need to choose our initial guess carefully
# should be between 2 and 3 by the "eyeball norm"
xk = 2.5 # initial guess - looks pretty close
tol = 1e-15 #tolerance to decide when to terminate
for i in range(100): #cap the number of iterations at 100
h = -rprime(xk)/r2prime(xk)# define step size
xk += h
if (abs(h) < tol): # if step size is small, this is our last iteration
print(i)
break
print(xk)
print(r(xk))
print(rprime(xk))
# In[ ]: