import numpy as np
# from enthought.mayavi import mlab
'''
ogrid[-1:5:6j,-1:5:6j]
[array([[-1. ],
[ 0.2],
[ 1.4],
[ 2.6],
[ 3.8],
[ 5. ]]), array([[-1. , 0.2, 1.4, 2.6, 3.8, 5. ]])]
'''
x,y = np.ogrid[-2:2:20j,-2:2:20j] #返回两个数组，一个长度为1，一个列数为1。前三个参数同range，但两边都是闭区间，当第三个参数是复数(j)时，以整数部分作步长。
z=x*np.exp(x**2-y**2)
#p1=mlab.surf(x,y,z,warp_scale="auto")
#mlab.axes(xlable='x',ylable='y',zlable='z')
#mlab.outline(p1)
#mlab.show()
#这个mayavi画图软件 太难安装了，遇到了错误，貌似需要安装对应的visual c然后设置visual c 的环境什么之类的才能好，暂时不管了
from scipy.optimize import leastsq
import pylab as p1
#最小二乘拟合(Least-Square fitting) ，已知数据集（x,y），已知函数y=f(x),找函数f中用的参数集p，使s(p)=∑[yi-f(xi,p)]?
def func(x,p):
'''数据拟合所用的函数：A*sin(2*pi*k*x + theta)'''
A,k,theta = p
return A * np.sin(2 * np.pi * k * x + theta)
def residuals(p,y,x):
'''实验数据x,y 和 拟合函数之间的差，p为拟合要找到的系数'''
return y - func(x,p)
x=np.linspace(0,-2*np.pi,100) #返回以0开始，以-2*np.pi结束的100个元素的等差数列
A,k,theta = 10,0.34,np.pi/6 #真实数据的函数参数
y0=func(x,[A,k,theta]) #真是数据
y1=y0+2*np.random.randn(len(x)) #加入噪声之后的实验数据
p0=[7,0.2,0] #第一次猜测的函数拟合参数
#调用leastsq进行数据拟合
#residuals为计算误差的函数
#p0为拟合参数的初始值
#args为需要拟合的试验数据
plsq=leastsq(residuals,p0,args=(y1,x))
print(u"真是参数：",[A,k,theta])
print(u"拟合参数：",plsq[0]) #实验数据拟合后的参数
p1.plot(x,y0,label=u"真实数据")
p1.plot(x,y1,label=u"带噪声的实验数据")
p1.plot(x,func(x,plsq[0]),label=u"拟合数据j")
p1.legend()
p1.show()
#对于一个离散的线性时不变系统h,如果它的输入是x,那么其输出y可以用x和h的卷积表示 :y=x * h
#现在的问题是如果已知系统的输入x和输出y，如何计算系统的传递函数h，或如果已知系统的传递函数h和系统的输出y，如何计算系统的输入x.这种运算被称为反卷积
#运算，是十分困难的，特别是在实际的运用中，测量系统的输出总是存在误差。
#本程序用各种fmin函数求卷积的逆运算
import scipy.optimize as opt
import numpy as np
def test_fmin_convolve(fminfunc,x,h,y,yn,x0):
''' x (*) h =y ,(*) 表示卷积
yn 为在y的基础上添加一些干扰噪声的结果
x0为求解x的初始值
'''
def convolve_func(h):
''' 计算yn - x (*) h 的 power
fmin 将通过计算使得此power最小'''
return np.sum((yn - np.convolve(x,h))**2)
#调用fmin函数，以x0为初始值
h0=fminfunc(convolve_func,x0)
print(fminfunc.__name__)
print('------------------')
#输出x (*) h0 和y 之间的相对误差
print('error of y:',np.sum((np.convolve(x,h0)-y)**2)/np.sum(y**2))
#输出h0和h之间的相对误差
print('error of h:',np.sum((h0-h)**2)/np.sum(h**2))
def test_n(m,n,nscale):
''' 随机产生 x,h,y,yn,x0 等数列，调用各种fmin函数求解
m 为 x 的长度，n 为 h 的长度，nscale 为干扰的强度 '''
x=np.random.rand(m)
h=np.random.rand(n)
y=np.convolve(x,h)
yn=y+np.random.rand(len(y))* nscale
x0=np.random.rand(n)
test_fmin_convolve(opt.fmin,x,h,y,yn,x0)
test_fmin_convolve(opt.fmin_powell,x,h,y,yn,x0)
test_fmin_convolve(opt.fmin_cg,x,h,y,yn,x0)
test_fmin_convolve(opt.fmin_bfgs,x,h,y,yn,x0)
if __name__=='__main__':
test_n(200,20,0.1)
#卡方检验示例代码
from scipy.stats import chi2
import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(1, 1)
df = 20
mean, var, skew, kurt = chi2.stats(df, moments='mvsk')
#绘制函数的起始点和终止点
#pdf为概率密度函数 chi2.pdf(x, df) = 1 / (2*gamma(df/2)) * (x/2)**(df/2-1) * exp(-x/2)
#默认的是标准（正态分布）概率密度函数，要设置它的形状可以通过参数loc和scale。chi2.pdf(x, df, loc, scale)等价于
#百分比函数(PPF) :the inverse of the CDF. PPF 函数和连续分布函数CDF相逆，
#ppf(0.01, df)表示取值小于哪个点的概率之和等于0.01,在df取值df,的情况下。
initial=chi2.ppf(0.01, df)
end=chi2.ppf(0.99, df)
print(chi2.cdf(end,df),chi2.cdf(initial,df))
#0.99 0.01 chi2.cdf(x,df)是chi2.ppf(x,df)的反函数。
#cdf (cumulative distribution function) 累计概率分布函数： 对于连续函数，所有小于等于a的值，其出现概率的和。F(a)=P(x<=a)
#随机变量小于等于某个数值的概率。F(x)为增函数。 累计分布的反函数，可以用来生成服从该随机分布的随机变量。
x = np.linspace(initial,end, 100)
#概率密度函数用于绘图
ax.plot(x, chi2.pdf(x, df), 'r-', lw=5, alpha=0.6, label='chi2 pdf') #‘r-’表示红色的实线，lw=5表示线宽度值取5，alpha=0.6表示透明度取0.6，label设置线的标注
ax.plot(x, chi2.pdf(x, 2), 'g-', lw=3, alpha=0.6, label='chi2 pdf')
ax.plot(x, chi2.pdf(x, 4), 'b--', lw=3, alpha=0.6, label='chi2 pdf')
ax.plot(x, chi2.pdf(x, 11), 'y--', lw=3, alpha=0.6, label='chi2 pdf')
#ax.axis([-3,2,-0.01,10]) #设置坐标轴的取值范围为x[-3,2],y[-0.01,10]
plt.title("df is %d"%df) #设置图形的标题
plt.show() #展示图形
rv = chi2(df) #绑定卡方分布的df值
plt.plot(x, rv.pdf(x), 'ko', lw=2, label='frozen pdf')
#plt.legend() #legend()用来设置图形标注框在图像中的位置
plt.show()
print(chi2.ppf([0.001, 0.5, 0.999], df))
print(chi2.ppf(0.001, df),chi2.ppf( 0.5, df),chi2.ppf( 0.999, df),sep=' @')
#[ 5.92104075 19.33742923 45.31474662]
#5.92104074549 @19.3374292294 @45.3147466181
r = chi2.rvs(df, size=100) #产生随机数100个
print(r)
#卡方公式(o-e)^2 / e
#期望值和收集到数据不能低于5，o(observed)观察到的数据，e（expected）表示期望的数据
#(o-e)平方，最后除以期望的数据e
import numpy as np
from scipy import stats
from scipy.stats import chisquare
list_observe=[30,14,34,45,57,20]
list_expect=[20,20,30,40,60,30]
print(chisquare(f_obs=list_observe, f_exp=list_expect))
p=chisquare(f_obs=list_observe, f_exp=list_expect)[1]
'''
返回NAN，无穷小
'''
if p>0.05 or p=="nan":
print("H0 win,there is no difference")
else:
print("H1 win,there is difference")
#Power_divergenceResult(statistic=11.441666666666666, pvalue=0.043293130315804972)
#H1 win,there is difference