Scribd
Upload
58 views15 pages

Simulation Code

This simulation code is performing the following: 1. Generating random samples and calculating the maximum likelihood (ML) and maximum a posteriori (MAP) estimators for different values of frequency and signal to noise ratio. 2. Computing the inverse variance of the estimators and comparing the performance of the ML estimators for different frequencies. 3. Comparing the performance of the MAP and ML estimators using the computed inverse variances.

Uploaded by

Sairam Gamidi
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
58 views15 pages

Simulation Code

This simulation code is performing the following: 1. Generating random samples and calculating the maximum likelihood (ML) and maximum a posteriori (MAP) estimators for different values of frequency and signal to noise ratio. 2. Computing the inverse variance of the estimators and comparing the performance of the ML estimators for different frequencies. 3. Comparing the performance of the MAP and ML estimators using the computed inverse variances.

Uploaded by

Sairam Gamidi
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
You are on page 1/ 15

SIMULATION CODE:

(a)Performance of ML estimators for diferent values of frequency ω0=0.7 and


ω0=0.05 with values of ϴ0=0.25π and N=11.

import csv
from matplotlib.pyplot import
show,hist,subplot,figure,plot,xlabel,ylabel,legend,title,xticks,yticks
from numpy import
sqrt,log,sin,cos,pi,random,linspace,arange,append,array,divide,transpose,matmul,li
nalg,dot,zeros,math,angle,reciprocal,log10
import matplotlib as plt
import numpy as np
from sklearn.metrics import r2_score
import math
import mpmath

a=10
w=0.05
pi=3.14
tita=pi/4
alpha=50
bita=50
n=11

def gauss(u1,u2):
z1=(sqrt(-2*log(u1))*cos(2*pi*u2))
z2=(sqrt(-2*log(u1))*sin(2*pi*u2))
return z1,z2

u1=random.uniform(0,1,n)
u2=random.uniform(0,1,n)

z1,z2= gauss(u1,u2)
#l=complex(z1,z2)
#print(z1)

#print(z2)

mag=[]
for i in range(len(z1)):
mag.append(abs(complex(a+z1[i],z2[i])))
#print(mag)
mag1=math.fsum(mag)

sum=0
angle1=[]
for k in range(n-1):
angle1.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z2[k]/mag[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])
sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z2[k]
sum1=sum1+z2[k]
sum2=sum2+k*mag[k]
sum3=sum3+k*k*mag[k]
p=sum
q=sum1
r=sum2
s=sum3

e=float(((p*mag1)-(q*r))/((s*mag1)-(r*r)))
t=(e*e)
#print(t)
m=1/(t)
print((m))

z3=z1*sqrt(100)
z4=z2*sqrt(100)
mag2=[]
for i in range(len(z1)):
mag2.append(abs(complex(a+z3[i],z4[i])))
#print(mag)
mag3=math.fsum(mag2)

sum=0
angle2=[]
for k in range(n-1):
angle2.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z4[k]/mag2[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z4[k]
sum1=sum1+z4[k]
sum2=sum2+k*mag2[k]
sum3=sum3+k*k*mag2[k]
p1=sum
q1=sum1
r1=sum2
s1=sum3

e1=float(((p1*mag3)-(q1*r1))/((s1*mag3)-(r1*r1)))
t1=(e1*e1)
#print(t)
m1=1/(t1)
print((m1))
z5=z1*sqrt(10)
z6=z2*sqrt(10)
mag4=[]
for i in range(len(z1)):
mag4.append(abs(complex(a+z5[i],z6[i])))
#print(mag)
mag5=math.fsum(mag4)

sum=0
angle3=[]
for k in range(n-1):
angle3.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z6[k]/mag4[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z6[k]
sum1=sum1+z6[k]
sum2=sum2+k*mag4[k]
sum3=sum3+k*k*mag4[k]
p2=sum
q2=sum1
r2=sum2
s2=sum3

e2=float(((p2*mag5)-(q2*r2))/((s2*mag5)-(r2*r2)))
t2=(e2*e2)
#print(t)
m2=1/(t2)
print((m2))

z7=z1*sqrt(20)
z8=z2*sqrt(20)
mag6=[]
for i in range(len(z1)):
mag6.append(abs(complex(a+z7[i],z8[i])))
#print(mag)
mag7=math.fsum(mag6)

sum=0
angle4=[]
for k in range(n-1):
angle4.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z8[k]/mag6[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z8[k]
sum1=sum1+z8[k]
sum2=sum2+k*mag6[k]
sum3=sum3+k*k*mag6[k]
p3=sum
q3=sum1
r3=sum2
s3=sum3

e3=float(((p3*mag7)-(q3*r3))/((s3*mag7)-(r3*r3)))
t3=(e3*e3)
#print(t)
m3=1/(t3)
print((m3))

sig=[50,100,10,20]
snr=[]
for i in sig:
snr.append((a*a)/i)
snr.sort()
#print(snr)
snrdb=[]

for i in snr:
snrdb.append(10*log10(i))

snrdb.sort()

print(snrdb)

a=10
w1=0.7
pi=3.14
tita=pi/4
alpha=50
bita=50
n=11

ma=[]
for i in range(len(z1)):
ma.append(abs(complex(a+z1[i],z2[i])))
#print(mag)
ma1=math.fsum(ma)

sum=0
angl1=[]
for k in range(n-1):
angl1.append(angle(complex(cos(w1*k+tita),sin(w1*k+tita)))
+math.asin(z2[k]/ma[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z2[k]
sum1=sum1+z2[k]
sum2=sum2+k*ma[k]
sum3=sum3+k*k*ma[k]
p4=sum
q4=sum1
r4=sum2
s4=sum3

e4=float(((p4*ma1)-(q4*r4))/((s4*ma1)-(r4*r4)))
t4=(e4*e4)
#print(t)
m4=1/(t4)
print((m4))

z3=z1*sqrt(100)
z4=z2*sqrt(100)
ma2=[]
for i in range(len(z1)):
ma2.append(abs(complex(a+z3[i],z4[i])))
#print(mag)
ma3=math.fsum(ma2)

sum=0
angl2=[]
for k in range(n-1):
angl2.append(angle(complex(cos(w1*k+tita),sin(w1*k+tita)))
+math.asin(z4[k]/ma2[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z4[k]
sum1=sum1+z4[k]
sum2=sum2+k*ma2[k]
sum3=sum3+k*k*ma2[k]
p5=sum
q5=sum1
r5=sum2
s5=sum3

e5=float(((p5*ma3)-(q5*r5))/((s5*ma3)-(r5)))
t5=(e5*e5)
#print(t)
m5=1/(t5)
print((m5))

z5=z1*sqrt(10)
z6=z2*sqrt(10)
ma4=[]
for i in range(len(z1)):
ma4.append(abs(complex(a+z5[i],z6[i])))
#print(mag)
ma5=math.fsum(ma4)

sum=0
angl3=[]
for k in range(n-1):
angl3.append(angle(complex(cos(w1*k+tita),sin(w1*k+tita)))
+math.asin(z6[k]/ma4[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z6[k]
sum1=sum1+z6[k]
sum2=sum2+k*ma4[k]
sum3=sum3+k*k*ma4[k]
p6=sum
q6=sum1
r6=sum2
s6=sum3

e6=float(((p6*ma5)-(q6*r6))/((s6*ma5)-(r6*r6)))
t6=(e6*e6)
#print(t)
m6=1/(t6)
print((m6))

z7=z1*sqrt(20)
z8=z2*sqrt(20)
ma6=[]
for i in range(len(z1)):
ma6.append(abs(complex(a+z7[i],z8[i])))
#print(mag)
ma7=math.fsum(ma6)

sum=0
angl4=[]
for k in range(n-1):
angl4.append(angle(complex(cos(w1*k+tita),sin(w1*k+tita)))
+math.asin(z8[k]/ma6[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z8[k]
sum1=sum1+z8[k]
sum2=sum2+k*ma6[k]
sum3=sum3+k*k*ma6[k]
p7=sum
q7=sum1
r7=sum2
s7=sum3

e7=float(((p7*ma7)-(q7*r7))/((s7*ma7)-(r7*r7)))
t7=(e7*e7)
#print(t)
m7=1/(t7)
print((m7))

sig=[50,100,10,20]
snr1=[]
for i in sig:
snr1.append((a*a)/i)
snr1.sort()
#print(snr)
snrdb1=[]
y=[]
for i in snr1:
snrdb1.append(10*log10(i))

snrdb1.sort()

print(snrdb1)

y=[]
m=[m2,m1,m,m3]
for i in m:
y.append(10*log10(i))
y.sort()

z=[]
m8=[m4,m5,m6,m7]
for i in m8:
z.append(10*log10(i))
z.sort()
#for i in range(len(y)):

print(z)
print(y)

plot(snrdb,y,marker='o',label='ML estimator for w=0.7')


plot(snrdb1,z,marker='*',label='ML estimator for w=0.05')
xticks(arange(0,10,2))
yticks(arange(15,60,5))
xlabel('snr in db')
ylabel('Inverse variance and ICRLB(db)')
title('performance comparision of ML estimators')
legend(loc='upper right')
show()

(b)Performance comparision between MAP and ML estimators with N=11 samples


import csv
from matplotlib.pyplot import
show,hist,subplot,figure,plot,xlabel,ylabel,legend,title,xticks,yticks
from numpy import
sqrt,log,sin,cos,pi,random,linspace,arange,append,array,divide,transpose,matmul,li
nalg,dot,zeros,math,angle,reciprocal,log10
import matplotlib as plt
import numpy as np
from sklearn.metrics import r2_score
import math
import mpmath

a=10
w=0.5
pi=3.14
tita=pi/4
alpha=50
bita=50
n=11

def gauss(u1,u2):
z1=(sqrt(-2*log(u1))*cos(2*pi*u2))
z2=(sqrt(-2*log(u1))*sin(2*pi*u2))
return z1,z2

u1=random.uniform(0,1,n)
u2=random.uniform(0,1,n)

z1,z2= gauss(u1,u2)
#l=complex(z1,z2)
#print(z1)

#print(z2)

mag=[]
for i in range(len(z1)):
mag.append(abs(complex(a+z1[i],z2[i])))
#print(mag)
mag1=math.fsum(mag)

sum=0
angle1=[]
for k in range(n-1):
angle1.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z2[k]/mag[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z2[k]
sum1=sum1+z2[k]
sum2=sum2+k*mag[k]
sum3=sum3+k*k*mag[k]
p=sum
q=sum1
r=sum2
s=sum3

e=float(((p*mag1)-(q*r))/((s*mag1)-(r*r)))
t=(e*e)
#print(t)
m=1/(t)
print((m))

z3=z1*sqrt(100)
z4=z2*sqrt(100)
mag2=[]
for i in range(len(z1)):
mag2.append(abs(complex(a+z3[i],z4[i])))
#print(mag)
mag3=math.fsum(mag2)

sum=0
angle2=[]
for k in range(n-1):
angle2.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z4[k]/mag2[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z4[k]
sum1=sum1+z4[k]
sum2=sum2+k*mag2[k]
sum3=sum3+k*k*mag2[k]
p1=sum
q1=sum1
r1=sum2
s1=sum3

e1=float(((p1*mag3)-(q1*r1))/((s1*mag3)-(r1*r1)))
t1=(e1*e1)
#print(t)
m1=1/(t1)
print((m1))

z5=z1*sqrt(10)
z6=z2*sqrt(10)
mag4=[]
for i in range(len(z1)):
mag4.append(abs(complex(a+z5[i],z6[i])))
#print(mag)
mag5=math.fsum(mag4)

sum=0
angle3=[]
for k in range(n-1):
angle3.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z6[k]/mag4[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z6[k]
sum1=sum1+z6[k]
sum2=sum2+k*mag4[k]
sum3=sum3+k*k*mag4[k]
p2=sum
q2=sum1
r2=sum2
s2=sum3

e2=float(((p2*mag5)-(q2*r2))/((s2*mag5)-(r2*r2)))
t2=(e2*e2)
#print(t)
m2=1/(t2)
print((m2))

z7=z1*sqrt(20)
z8=z2*sqrt(20)
mag6=[]
for i in range(len(z1)):
mag6.append(abs(complex(a+z7[i],z8[i])))
#print(mag)
mag7=math.fsum(mag6)

sum=0
angle4=[]
for k in range(n-1):
angle4.append(angle(complex(cos(w*k+tita),sin(w*k+tita)))
+math.asin(z8[k]/mag6[k]))
#t = angle(angle1, deg=True)
#print(angle1[k])

sum=0
sum1=0
sum2=0
sum3=0
for k in range(n-1):
sum=sum+k*z8[k]
sum1=sum1+z8[k]
sum2=sum2+k*mag6[k]
sum3=sum3+k*k*mag6[k]
p3=sum
q3=sum1
r3=sum2
s3=sum3

e3=float(((p3*mag7)-(q3*r3))/((s3*mag7)-(r3*r3)))
t3=(e3*e3)
#print(t)
m3=1/(t3)
print((m3))

sig=[50,100,10,20]
snr=[]
for i in sig:
snr.append((a*a)/i)
snr.sort()
#print(snr)
snrdb=[]
y=[]
for i in snr:
snrdb.append(10*log10(i))
y.append(10*log10(i))
snrdb.sort()

print(snrdb)
snr1=[1,3,7,10]
y=[]
m=[m2,m1,m,m3]
for i in m:
y.append(10*log10(i))
y.sort()

f=mpmath.besselj(0,50)
#print(f)
sum=0

for i in range(1,5000):
sum=sum+((((-1)**i)*mpmath.besselj(i,50))/(i**2))
#print(sum)

sig1=(pi*pi)/3+(4/f)*(sum)
c=(n*(n-1)*((2*n)-1))/6
d=(n*n)*((n*n)-1)/3
#print(sig1)
#print(d)
#print(c)
P=[]
for i in snr1:
P.append(2*i*(c*bita*bita+i*d*bita+(i*i*d)))
#print(P)
pos=[]
p=[]
s=[]
for i in range(len(snr1)):
p.append(((((2*snr1[i]*alpha*n+alpha*bita)**2)*(sig1))+(((snr1[i]*bita*n*(n-
1))**2)*(sig1))+P[i]))
s.append(((snr1[i]*2*c+alpha)*(2*snr1[i]*n+bita)-(snr1[i]*n*(n-1))**2)**2)
#pos.append(((((2*i*alpha*n+alpha*bita)**2)*(sig1))+(((i*bita*n*(n-
1))**2)*(sig1))+P[i])/((i*2*c+alpha)*(2*i*n+bita)-(i*n*(n-1))**2)**2)
pos.append(float(p[i]/s[i]))
print(pos)
#print(s)
#print(p)
pos1=[]
for i in pos:
pos1.append(1/i)
print(pos1)
z=[]
for i in pos1:
z.append(10*log10(i))
print(z)
print(y)
plot(snrdb,y,marker='o',label='ML estimator')
plot(snrdb,z,marker='*',label='MAP estimator')
xticks(arange(0,10,2))
yticks(arange(15,60,5))
xlabel('snr in db')
ylabel('inverse variance/mse(db)')
title('performance comparision of MAP and ML estimators')
legend(loc='upper right')
show()

SIMULATION RESULTS:
Fig. ML estimator performance for values of frequency ω0=0.7 and ω0=0.05
with values of ϴ0=0.25π and N=11.
Fig. Performance comparision between MAP and ML estimator with N=11
samples

CONCLUSION:

The comparision shows MAP estimator has a better performance gain


over ML estimator due to incorporating apriori knowledge about the
parameters.Also as SNR increases the difference of gain between MAP and ML
estimators decreases.The performance of ML estimators for diferent values of
frequency ω0=0.7 and ω0=0.05 is observed to be approximately equal at high
SNR for values of ϴ0=0.25π and N=11.

You might also like